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PIECE-WISE  ANALYTIC  EVALUATION  OF  THE  RADIATIVE  TAIL 
FROM  ELASTIC  AND  INELASTIC  ELECTRON  SCATTERING 
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National  Bureau  of  Standards 
Washington,  D.C.  20234 

S.  E.  Williamson 

Nuclear  Physics  Laboratory,  Department  of  Physics 
University  of  Illinois  at  Urbana-Champaign 
Urbana,  II 1 1 nois  61801 

ABSTRACT 

We  consider  here  the  calculation  of  the  radiative  tail  from  the  elastic 
peak  in  medium  and  high  energy  electron  scattering  as  well  as  from  a discrete 
inelastic  level  of  the  recoiling  nucleus.  We  examine  the  method  generally 
used  for  this  calculation,  viz.,  a numerical  integration  of  the  differential 
cross  section  over  the  angles  of  the  unobserved  photon,  and  discuss  the 
difficulties  inherent  in  this  numerical  integration  due  to  the  sharp  peaking 
of  the  integrand.  We  present  an  alternative  method  for  calculating  the 
radiative  tail,  in  which  the  region  of  integration  is  divided  into  an 
arbitrary  number  of  subintervals,  the  structure  functions  are  fitted  by  cubic 
spline  functions  in  each  subinterval,  and  the  integrations  are  then  performed 
analytically  in  closed  form.  This  method  has  the  advantages  of  greatly 
increased  accuracy  and  a reduction  of  the  computation  time  by  a factor  which 
can  vary  between  10  and  103,  depending  on  the  kinematics. 


Key  words:  electron  scattering  cross  sections;  elastic  electron  scattering; 

elastic  radiative  tail;  inelastic  electron  scattering;  inelastic 
radiative  tail;  intermediate  and  high  energy  electron  scattering. 
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I.  Introduction 


The  analysis  of  medium  and  high  energy  electron  scattering  experiments 
involves,  as  is  well  known,  the  calculation  of  electromagnetic  processes 
(e.g.,  radiative  corrections,  Landau  straggling,  radiative  tails,  etc.)  which 
inevitably  accompany  and  mask  in  some  manner  the  electron-nucleus  scattering 
of  primary  interest.  The  methods  for  evaluating  the  cross  sections  for  these 
processes  has  been  given  in  numerous  articles  and  reviews  familiar  to 
experimenters,  and  are  generally  incorporated  into  the  analysis  of  electron 
scattering  experiments.  We  consider  here  the  calculation  of  the  radiative 
tail  from  the  elastic  peak,  as  well  as  from  a discrete  inelastic  level  of  the 
recoiling  nucleus.1  We  examine  the  method  which  is  generally  used  for  this 
calculation,  namely  a numerical  integration  of  the  differential  cross  section 
over  the  angles  of  the  unobserved  photon,  and  discuss  the  difficulties 
inherent  in  this  numerical  integration.  Finally,  we  present  an  alternative 
method  for  calculating  the  radiative  tail  having  the  advantage  of  greatly 
increased  accuracy  and  a reduction  of  the  computation  time  by  a factor  which 
can  vary  between  10  and  103,  depending  on  the  kinematics  (see  Timing  Table, 
p.  56). 

The  formalism  used  for  the  calculation  of  the  radiative  tail  in  almost 
all  analyses  of  high-energy  electron  scattering  experiments  is  that  presented 
in  a detailed  series  of  papers  by  Tsai  [1-3] 2 and  later  refined  by  Miller  [4], 
In  all  of  these  papers  the  cross  section  for  the  scattering  of  an  electron  of 
initial  momentum  Ip^l,  final  momentum  |£2|  at  an  angle  6,  with  emission  of  a 
photon,  is  given  by  an  integral  of  the  form 


^he  assumption  that  the  final  state  of  the  recoiling  nucleus  is  a discrete 
level  does  not  constitute  a significant  restriction  in  practice,  since  in  the 
analysis  a final  continuum  state  is  generally  treated  as  a sum  of  many 
discrete  levels. 

2Figures  in  brackets  indicate  literature  references  at  the  end  of  this  paper. 


(1) 


dnd|£2l  = ^{w2f 2 + wif  i}  dnk  * 

Here  the  functions  f2  and  are  known  analytic  functions  of  the  initial  and 
final  electron  momenta  and  scattering  angle  (which  remain  fixed)  and  the 
photon  angles  0^,^,  over  which  one  integrates.  The  structure  functions  W2 
and  W1  are,  for  elastic  scattering,  a function  of  the  four -momentum  transfer 
q2  < 0.  They  are,  in  general,  qiven  as  the  result  of  a numerical  calculation 
in  which  a charge  distribution  has  been  entered  into  a phase  shift 
calculation.  (That  is,  although  a Born  approximation  formalism  is  used  to 
derive  the  cross  section  given  in  (1),  one  uses,  in  the  integrand  in  that 
expression,  "experimental"  structure  functions,  defined  either  by  the 
experimental  cross  section  divided  by  the  Mott  cross  section  or  by  a phase 
shift  calculation.)  It  is  expected,  therefore,  that  the  calculation  of  the 
radiative  tail  from  (1)  requires  a numerical  integration.  However,  the 
numerical  integration  of  the  expression  appearing  in  the  integrand  in  (1)  is 
not  completely  trivial.  Depending  on  the  kinematics,  the  functions  f2  and  f1 
can  be  very  sharply  peaked,  — a reflection  of  the  fact  that  the  photons  are 
emitted  preferentially  in  the  direction  of  either  the  initial  or  the  final 
electron.  In  addition,  since  the  functions  f2  and  fx  have  an  overall  factor 
— , the  integrand  may  rise  sharply  for  q2  near  its  minimum  value,  achieved  (in 

q4 

the  coordinate  system  defined  following  eqs  (9))  when  cos  0^  = 1.  These  peaks 
are  discussed  in  detail  both  in  Maximon  and  Isabelle  [5]  and  in  Mo  and  Tsai 
[2],  In  general,  the  height  of  the  peaks  is  of  order  (E/mc2)  relative  to  the 
background  and  the  width  is  of  order  (mc2/E)  relative  to  the  width  of  the 
region  over  which  one  integrates.  The  contribution  from  the  peaks  is  thus  of 
the  same  order  of  magnitude  as  that  from  the  background,  and  for  high  energies 
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great  care  must  be  taken  to  ensure  that  the  numerical  integration  procedure 
has  indeed  converged  and  gives  the  desired  accuracy.  The  result  of  the 
peaking  is  thus  that  the  numerical  integration  may  be  quite  lengthy,  or  of 
dubious  accuracy,  or  both.  Further,  given  the  large  number  of  intervals  that 
must  be  chosen  for  the  integration  in  the  region  of  the  peaks,  the  numerical 
evaluation  of  the  structure  functions  W2  and  W1  in  each  integration  interval 
would  be  prohibitive  in  computation  time,  and  rather  useless  since  the 
functions  W2  and  Wx  vary  little  over  even  the  entire  region  of  a peak.  The 
procedure  generally  followed3  is  therefore  to  divide  the  entire  region  of 
integration  into  a fixed  number  of  subintervals  (1  < i < N)  (say  with  N of  the 
order  of  100)  and  to  fit  W2  and  VJ1  by  spline  functions  in  each  of  these 
subintervals,  e.g.,  by  a cubic  function  of  -q2: 


W2  = ai  + bi  (-q2)  + c i (-q2)  2 + d^-q2)3 
Wi  = ai  + bi  (-q2)  + ci  (-q2)  2 + d^-q2)3  . 


(2) 


With  the  overall  region  of  integration  divided  into  these  subintervals  (some 
of  which  contain  the  peaks)  and  W2  and  WL  represented  by  these  simple  analytic 
expressions  in  each  subinterval,  the  numerical  integration  is  then  carried  out 
as  just  discussed,  the  peak  rendering  the  integration  costly  in  time  and 
questionable  in  accuracy. 


3This  is,  essentially,  the  procedure  followed  by  R.  Altemus  in  the  analysis  of 
the  data  on  56Fe  taken  at  Bates  (reported  in  Phys.  Rev.  Lett.  44,  965 
(1980)),  by  M.  Deady  in  the  analysis  of  the  data  on  40,  48Ca  taken  at  Bates 
(reported  in  Phys.  Rev.  C28,  631  (1983)),  and  by  Z.  Meziani  in  the  analysis 
of  the  data  on  40’  48Ca  and  56Fe  taken  at  Sacla.y  (to  be  published  shortly). 
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We  note  that  the  integration  procedure  just  described  consists  in  fact  of 
two  parts  which  are  quite  independent.  In  the  first,  each  of  the  structure 
functions  is  fitted  by  a spline  (or  other  useful  function)  in  each  of  the  N 
subintervals  into  which  the  integration  region  is  subdivided  for  this  purpose. 
(For  example,  for  a point  charge  with  spin  zero  we  have  W2  = Z2  and  W1  = 0,  so 
that  we  may  choose  N = 1.)  In  the  second  part,  the  numerical  integration  of 

(1)  is  performed  in  each  of  the  N subintervals,  with  W2  and  WL  now  given  by 

(2) ,  using  some  suitable  integration  technique  (e.g.,  Simpson's  rule  or 
Gaussian  integration).  In  the  event  that  a peak  occurs  in  a given  sub- 
interval, then  clearly  this  subinterval  must  be  further  subdivided,  and  the 
convergence  of  the  integration  procedure  must  be  tested. 

The  essence  of  the  present  report  is  contained  in  the  observation  that 
once  we  have  carried  out  the  fitting  of  W2  and  WL  by  functions  of  the  form 
given  in  (2)  (the  first  part  of  the  procedure  just  described),  there  is  no 
longer  any  need  to  perform  the  integration  numerically.  With  W 2 and  Wx  given 
by  finite  polynomials  in  q2  the  functions  f2  and  f1  are  sufficiently  simple 
that  the  integration  can  be  performed  in  each  subinterval  in  closed  form. 

(The  resulting  functions  involve  only  algebraic  functions  and  logarithms.) 

The  integration  over  the  entire  region  then  consists  merely  of  summing  the 
contribution  from  each  of  the  N subintervals.  The  problems  related  to  these 
peaks  are  thereby  avoided  completely.  The  presence  of  sharp  peaks  does  not 
necessitate  a further  division  of  the  subintervals.  It  is  thus  by  eliminating 
the  numerical  integration,  the  second  part  of  the  procedure  usually  followed, 
that  we  reduce  the  time  for  the  entire  calculation  by  a factor  which  varies 
between  10  and  103,  depending  on  the  specific  kinematics.  The  accuracy  of  the 
entire  calculation  is  also  greatly  increased,  since  it  is  now  given  by  the 
precision  of  the  computer  rather  than  the  accuracy  of  the  particular  numerical 
integration  method. 
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As  a last  general  comment  before  entering  into  the  details,  we  observe 
that  in  this  note  the  method,  which  we  call  "piece-wise  analytic"  evaluation, 
is  applied  to  the  integration  of  the  doubly  differential  cross  section  over 
the  angles  of  the  photon.  However,  it  is  equally  applicable  (with,  of  course, 
different  resulting  expressions)  to  the  situation  of  interest  for  the  analysis 
of  experiments  using  photons  as  the  electromagnetic  probe.  Then  one  inte- 
grates this  same  doubly  differential  cross  section  over  the  angles  of  the 
electron:  The  expressions  which  result  in  this  latter  case  will  be  presented 

in  a future  report,  in  which  we  consider  in  particular  the  cross  section  for 
polarized  photons  (integrated  over  the  angles  of  the  final  electron). 

We  next  present  the  details  of  the  calculation  of  the  radiative  tail 
given  in  this  report.  We  divide  this  presentation  into  three  parts:  In  the 

first  we  give  preliminaries  to  the  method  presented  here.  Notation  is 
defined,  the  differential  cross  section  which  forms  the  basis  of  our  work  is 
presented,  and  a few  comments  are  made.  In  the  second  we  give  the  details  of 
the  calculation  as  reflected  directly  in  the  computer  program,  a copy  of  which 
is  included  with  this  report.4  Finally,  in  the  Appendix  we  give  details  of 
the  algebra  involved  in  obtaining  our  final  expression  for  the  cross  section, 
given  by  eq  (38)  and  the  subsidiary  equations  defining  the  quantities  in  it. 


4The  Fortran  program  is  available  from  the  authors.  Requests  should  be 
accompanied  by  a magnetic  tape  and  a specification  of  the  desired  format. 
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II.  Preliminaries  for  the  Piece-Wise  Analytical  Evaluation  of  the  Radiative 


Tail 


The  cross  section  under  consideration  here  is  the  radiative  tail  to  both 
elastic  and  inelastic  scattering  of  electrons  from  a nucleus.  In  this  report 
(as  in  the  work  of  Tsai)  we  consider  emission  of  photons  only  by  the  electron. 
Furthermore,  the  formalism  used  is  that  of  the  Born  approximation.  Thus  in 
terms  of  diagrams  we  have 


We  do  not  consider  here  the  contribution  due  to  radiation  by  the  recoiling 
nucleus.  That  has  been  treated  in  an  approximate  manner  by  Miller  [4]  and 
will  be  examined  by  us  in  a report  to  appear  later. 

The  notation  used  here  is  as  follows:  Energies,  momenta,  and  masses  are 
all  understood  to  be  in  energy  units.  When  four-vectors  are  written,  the 
metric  is  as  in  Bjorken  and  Drel 1 [6],  viz., 


V 


Y 


(3) 
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Pi  = ( £1\P  l) : 

P2  = (e2>P2): 
k = ( w,J< ) : 

q = (q0»^): 

P!  ■ l.il)  = 

* (M,0) 

P2  = (E2’£.2>: 

= (M  + q0,q) 

m = 

M = 

* 

M = 


four-momentum  of  the  incident  electron 

four -moment urn  of  the  scattered  electron 

four-momentum  of  the  emitted  (real)  photon 

four -momentum  of  the  virtual  photon  exchanged  between 

the  electron  and  the  target  particle 

four -momentum  of  the  target  particle 

in  the  lab  system 

four -momentum  of  the  recoiling  nucleus 

in  the  lab  system 

rest  mass  of  the  electron 

rest  mass  of  the  target  particle. 

rest  mass  of  the  final,  recoiling  nucleus 


(=  (M*  - M)  Cl  + *4^)  ” M*  - M for  144*  « 1).  (4) 


Then  from  energy-momentum  conservation 


pi  + Pi  = p2  + P2  + k • 


(5) 


The  momentum  transferred  to  the  nucleus  is  then 


q = P2-P1  = P1-P2-k  . 


(6) 


Further,  it  is  useful  to  define  the  four-vector  u by 


u = (u  0,_u)  = p1  + P1-p2  = P2+k 


(7) 
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Then  in  the  lab  system 


U0  = M + Ej  - c2 


(8a) 


and 


£ = £1  - £2 


(8b) 


We  note  that  u0  and  £ are  independent  of  both  photon  energy  and  angles. 

0 = scattering  angle  of  the  electron  (0  = $ (jP  1*P.2) ) 

0^  = angle  between  £ and  k_ 

0^  = angle  between  £ and  £x 


As  in  Mo  and  Tsai  [2]  the  integration  over  photon  angles  is  performed  in  a 
coordinate  system  such  that  the  z axis  is  in  the  direction  of  £ = £:  - £2.  As 
observed  there,  this  coordinate  system  has  the  particular  advantage  that  q2  is 
then  independent  of  the  photon  azimuthal  angle  . The  structure  functions 
W2(-q2)  and  W^-q2)  can  then  be  taken  outside  of  the  integral  over  which 
can  then  be  performed  analytically.  We  show  this  explicitely  in  the  following 
paragraph. 

From  the  notation  just  introduced,  since  the  recoiling  nucleus  is  assumed 


(9) 


*• 

to  be  a free  particle  of  mass  M , we  have 
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(10) 


M*2  = P22  = (u  - k)2 
= u2  - 2k  *u 


from  which 


0) 


- IjJ  | cos  6k  J ‘ 


(ID 


Thus  the  photon  energy,  w,  depends  on  the  photon  direction,  but  in  the 
particular  coordinate  system  we  have  chosen  it  is  a function  only  of  0^,  and 
not  of  (Note,  however,  from  the  expression  for  w given  below,  viz., 

Q 2 Q2 

w = - e£  - y + that  in  the  limit  ♦ 0,  we  have  w = ei  - e2  * v;  that 

is,  in  the  absence  of  recoil  the  photon  energy  is  independent  of  photon 
direction.) 

Similarly,  we  can  write 


M*2  = P22  = (pt  + q) 2 


+ q2  + 2q-Pj  . 


Thus  in  the  lab  system 


u . 


(12) 


(13) 
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On  the  other  hand,  substituting  q = - p2  - k in  2q*P1  above  gives 


(14) 


so  that  in  the  lab  system 


(15) 


Thus  in  the  coordinate  system  with  z-axis  in  the  direction  of  u.,  both  u>  and 
q2  are  independent  of  Only  the  products  p^k  and  p2*k  depend  on  <a  . 

Our  starting  point  is  then  the  doubly  differential  cross  section  for 
scattering  of  an  electron  with  the  emission  of  a photon,  essentially  as  qiven 
by  Tsai  [1]  and  Mo  and  Tsai  [2],  We  make  two  minor  changes  in  the  expressions 
presented  in  these  references,  discussed  immediately  following  the  cross 
section  written  below. 


d2o  a3  |^_2 1 2 f+1  2Mco  d (cos  0^ ) r 2ir 


dfld|£2|  (2")2  ejpj  J-1  q4(u0  - |_u  |cos  6k ) J° 


where 


I } * W2(-q2)  [-  [2«z(e  -»)  + %q2]  - [2e,(«2+»)  + iq2] 

L <pl'k>  (p2-k)2 

+ ipj-lofivkT  Cmz(p1*p2- “2)  + (Pi,P2)[2eiVPi‘P2t“(£r£2)]) 

- (p  'l-k-  [2(e1erc2We12)  + %q 2 - P^Pj  - m2] 

+ (p  -k)  [2(eiVel"fe22)  + *p2  - p i*p2  - m2] 

+ W,(-q2)  |m2(2m2+q2)  ( 1 + 1 ) 

L V ( p ! -k  ) 2 ( p 2*k ) 2/ 

4p  1*P2(p  1*p2-2m2) 

+ (p1*k)(p2*k) 

- (t^FT  * T^)(2pi’P2+2m2-q2) 

+ qj  (i7) 

The  first  modification  we  have  made  in  the  expression  of  Tsai  concerns  the 
structure  functions.  We  use  W2  and  Wx,  defined  by  the  elastic  scattering 
cross  section,  which  at  high  energies  is  given  in  terms  of  these  functions  by 

da  a2cos2^6[W9  + 2tan2^0  W,] 

— ^ = .--9- — . (18) 

dft  4e  2sin4^0  [l  + £_1  s i n 2% 0 ] 

M 


11 


(Note  that  for  a point  charge  with  spin  zero,  W2(0)  = Z2,  W^O)  = 0.)  This  is 

the  same  as  the  change  of  notation  introduced  by  Stein  et  al.  [7]  in 

Appendix  A of  their  paper.  That  appendix  is  a useful  reference  in  that  it  is 

essentially  a summary  of  Tsai's  SLAC-PUR  848  [3].  The  structure  functions  F. 

J 

and  G.  used  by  Tsai  are  related  to  W2  and  Wx  by 
j 

Fj  * 

G.  = 4M2W1  . (19) 

J 


The  second  modification  made  in  the  expression  of  Tsai  concerns  the  term 


ip2i2 

. , which  appears  as  a factor  to  the  integral  just  written.  In  the 

e2  l£.i  I 

articles  of  Tsai  this  factor  has  been  replaced  by  its  high  energy  limit,  viz., 
e2 

— . While  this  replacement  is  clearly  of  little  significance  numerically  for 

ei 

ip2i2 

energies  of  experimental  interest,  we  retain  the  factor  — r— -t  out  of  a sense 

e2 1J1  1 ' 

of  consistency,  since  no  high  energy  approximations  are  made  anywhere  else  in 
the  expression  for  the  cross  section.  Our  expression  is  then  identical  to 
that  given  by  Nguyen-Ngoc  and  Perez-y-Jorba  [8],  apart  from  a trivial  differ- 
ence in  the  normalization  of  the  structure  functions:  they  use  F and  G,  which 
are  related  to  our  W2  and  Wx  by 


F 


M2 


G = - W1  . 


(20) 
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In  any  event,  this  term  is  closely  related  to  the  physics:  Writing 

I P.2 1 2 /IP2I\  l£2l  IP2I 

— ; — r = — r,  we  note  that  the  factor  comes  from  the  fact  that  we 

£2Ip.ii  V e2  / Ipj  ; 

a 


have  written  the  cross  section  as 


d2 


j-  rather  than  . The  term  |pL 


dftd  I P2i  ^ “2 

comes  from  the  incident  flux,  and  the  final  factor  |p2|  comes  from  the  density 

of  final  states  for  the  scattered  electron:  (Note  that  in  Stein  et  al . [7] 

d2o  -2^ 

they  write  d^-e -,  so  that  only  should  appear,  but  they  replace  this  with 
e2 

— . ) This  factor  without  high  energy  approximation  has  also  been  retained  to 
ei 

facilitate  comparison  with  the  Bethe-Heitler  cross  section  in  the  limit  M 

Finally,  we  note  that  in  Stein  et  al.  [7]  the  structure  functions  W2  and 
Wx  have  been  replaced  by  W2  and  Wx,  defined  by  (A42)  and  (A43)  of  Stein 
et  al.: 


W_(  q 2 ) = F(q2)W,(q2) 


WJq2)  = F(q2)M1(q2)  (21) 

to  account  for  both  thick  target  bremsstrahlung  (equivalent  radiator)  and  that 
part  of  the  radiative  correction  which  does  not  depend  on  the  resolution  AE. 

In  addition,  in  the  work  of  Miller  [4]  these  factors  are  further  modified  to 
account  (in  an  approximate  manner)  for  the  multiple  emission  of  soft  photons 
from  both  the  electron  and  the  target  particle.  The  important  point  to  retain 
here  is  that  all  these  factors  depend  on  the  photon  direction  only  through  the 
four-momentum  transfer  q2  (or  via  the  photon  energy,  which  is  again  a linear 
function  of  q2,  as  we  have  just  seen).  All  these  effects  can  therefore  be 
accounted  for  simply  by  an  appropriate  modification  of  the  structure  functions 
which  leaves  them  independent  of  . 
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We  next  perform  the  integration  over  as  in  Mo  and  Tsai  [2],  Tsai  [3] 
and  Stein  et  al . [7],  and  obtain 


d2o  a3  |£2I2  f+1  2Ma)d(cos0u) 


/: 


dnd |p2|  2v  J~1  q4(u0  - McoseJ 


( 1 


(22) 


where 


f 1 * W2(-q2) 


am2  rn  _ , _ . \ . i a'm2 


[2£i(  e2+“)  + - , [2£2<£i-“)  + 


+ 2v(i  - y){m2(p1*p2-t*32)  + (p1*P2)[2e1e2-Pi,P2+aj(ei-e2)H 

+ eie2+eia>fe22)  + - P 1 *P 2 - m2] 

’ ^[2(  eie2-e2a>fei2)  + ^2  " P 1 *P  2 ' m2l 


-] 

+ W^-q2)  j^(j3  + p)  m2(2m2+q2)  + 4v  (■£  - -)(p  1*p2)  (p  1*p2-2m2) 
+ (j  - '“)(2p  1*p2+2m2-q2)  + ^ 


(23) 
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in  which5*6 


a = w(  e2  - |£2|cos0pcos0k) 
a'=  o)(  ex  - |£1|cos0scos0(<) 
b = - w|£2|sin0psin0k 
v = (a'  - a)-1 
x = (a2  - b2) 1/2 


y = (a1  2 - b2) 1/2 


|£x  | COS  0 - |£2| 

COs9p  ■ iuj 


cos  0„  = 
s 


l£xl  - I P.2 1 eos  0 


(24) 


5The  symbol  x which  is  defined  in  (24)  and  used  in  (23)  is  employed  at  this 
point  out  of  consistency  with  the  notation  of  Tsai  [3]  and  Stein  et  al . [7]. 
It  is  used  nowhere  else  in  this  report  as  defined  in  (24).  Elsewhere 
throughout  this  report  we  define  x = -q2. 

6The  expressions  given  here  in  (24)  are  those  given  in  Tsai's  SLAC-PUB  848  [3] 
and  Stein  et  al . [7],  They  represent  a considerable  simplification  of  the 
expressions  given  originally  in  Mo  and  Tsai  [2], 
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III.  Details  of  the  Piece-Wise  Analytical  Evaluation  of  the  Radiative  Tail 


At  this  point  we  depart  from  the  procedure  followed  in  these  references. 
Since  we  perform  the  remaining  integration  analytically,  the  most  appropriate 
variable  of  integration  is  q2  rather  than  cos  9^,  in  view  of  the  functions 
W2(-q2)  and  W^-q2).  Indeed,  the  form  of  the  integrand  itself  is  then 
somewhat  simpler.  To  give  just  one  example,  referring  to  the  overall  factor 
in  the  integrand  in  above  expressions  for  the  radiative  tail,  we  now  have, 
with  q2  as  integration  variable. 


We  will  not,  however,  encumber  this  section  of  the  report  with  these  algebraic 
manipulations.  They  are  presented  at  the  end  in  the  Appendix.  Here  we  give 
only  the  final  form  of  our  expression  for  the  radiative  tail,  which  lends 
itself  directly  to  the  computer  implementation  of  the  calculation. 

As  variable  of  integration  we  choose  finally  a positive  quantity,  viz.. 


2M  a)  d ( cos 0^)  1 d(q2) 


(25) 


q4(uQ  - |jj | cos 9k)  \u\  q4 


x 


-q2  > 0 . 


(26) 


The  cross  section  for  the  radiative  tail  is  then  given  by 


d2a  a3  |jd  2 1 2 

dftd|p_2|  tt  e^pj 


(27) 
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where 


1 1=  - [ 60  + SjX  + 62x 2] 


x 2D^  3/2(x) 


x 2D2  3 /2( x ) 


[«o  + 51X  + {2x2] 


|u|x2  (2A  - x)  VDj  l/2(x)  021/2(x), 


x2Dj  1/2(x)  x2D21/2(x) 


xD  1/2(x)  xD21/2(x) 


(28) 


and  where  j }'  is  obtained  from  { } merely  by  replacing  the  parameters 

30,31,...  by  the  corresponding  parameters  with  primes: 

{30,  S2,  60,  61#  52,  n,  Y,  Tp  t2,  pp  p2} 

♦ f V » V > e2' , 60‘ , «i' , 62‘  , n' , Y' , Ti' , t2'  , pi' , p2' } 


In  (28)  the  quadratic  functions  D^(x)  and  ( x)  are  related  to  the  quantities 

a,  a',  and  b of  (24)  by  eqs  (A19)  and  (A25).  Referring  to  the  terms  in  the 

1 


integrand  with  factor 


(2X  - x) 

indeed  lies  in  the  range  of  integration: 


, we  note  in  the  appendix  that  the  point  x = 2 A 


x < 2A  < xM  , 
m M 


(29) 
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but  that  D|(2A)  = (2 X) , so  that  there  is  in  fact  no  singularity  in  the 

integrand  at  the  point  x = 2A.  (This  is  shown  explicitely  later  in  this 
report,  in  the  Appendix.) 

The  parameters  appearing  in  the  integral  in  (27)  and  (28)  are  defined 
shortly  below.  However,  we  first  outline  our  piece-wise  analytic  evaluation 
of  the  cross  section.  We  divide  the  range  of  integration  (x  iX^)  into  N 
subintervals  (which  can  be  of  different  and  arbitrary  length)  writing 


(30) 


with 


x0  = xrn  - XN  = XM  ' (31) 

In  each  subinterval  (x-j_i»x-j)  we  represent  W2(x)  and  WL(x)  by  a third  order 
polynomial  in  x:  The  coefficients  of  the  polynomial  will,  in  general,  be 

different  in  each  subinterval. 

!n  (xi_1,xi) 


W2(x)  = a-  + b-x  + c.jX2  + d-x3 


I I 


Wx(x)  = a-  + b-x  + c-x2  + d-x3 


(32) 


We  may  then  write 
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/*xm  n rxi 

!j  dx[w2(x){  } + Wjx))  } ] = M dx[W2(x){  } + W x(x)  { } ] 

*m  i=1Vl 


N fXi  r 

vv*  dx[(ai 

1-1  xi-l 


bix 


+ C-X 


dix 


;)  { 


+ (a-  + brx  + c^x2  + d^x3)  { 


Considering  the  expressions  just  given  for  { } and  { } (the  functions  D^x) 

and  (x ) are  quadratic  functions  of  x),  the  integrals  can  now  all  be 
performed  analytically.  The  integrations  are  all  straightforward  — none 
result  in  functions  any  more  complicated  than  the  logarithm.  However,  care 
must  be  taken  in  arranging  the  terms,  and  in  particular  the  arguments  of  the 
logarithms,  so  that  at  high  energies  there  are  no  cancellations  between  almost 
equal  terms,  leaving  a remainder  of  relative  order  (m/ej)2,  with  a consequent 
loss  in  precision. 

Referring  to  the  integrand,  the  integrals  to  be  performed  analytically 
are  thus  of  the  the  form 


(1)  * P 

U Jv 


dx 


xi-l  x2^i3/2(x) 


[a^  + b^x  + c-x2  + d^ x 3 ] [e0  + $xx  + 62x 2 ] (34a) 


'i2>  = / - h + v ♦ CiX2  + d, 


x 3] 


(34b) 


xi-l  x 


is -r 


dx 


Xi.j  x2Dj1/2(x) 


[a.  + b^x  + c-x2  + d - x 3 ] 


(34c) 
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(4)  . f1 

■l.i  ' l 


dx 


xi-l  xDi1/2(x) 


[a-  + + ^x2  + d-x3] 


(34d) 


and 


■ f 


dx 


xi_l  (2X  - x)  \p11/2(x)  D2I/2(x)> 


[a-  + b-x  + c^2  + d^x3].  (34e) 


The  integrals  iP],  iP],  and  il4]  follow  from  iP],  iP],  and  l{4!  merely  by 

U.9I  L.  ^ I l j I J.  ^ I J.  ^ I 1 ) I 

replacing  D^(x)  by  D^(x)  in  the  integrands  (and  replacing  30,$1,62  in  I jP  by 
60 * 51 » 52 ^ • We  then  write 


I-  = 
1 


. t(1)  . t(1)  . _J]_  r(2)  + T r ( 3 ) 

!l,i  A2,i  |u|  M Vl.i 


x i(3) 
T212,  i 


P I(4)  + p I(4)  - yI(5) 
pi1l,i  ^2,1  Y1i 


(35) 


and  define  1^  by  substituting  throughout  the  corresponding  parameters  with 
primes: 

(30,  Bp  ... } + {Bq,B^,...} 
and 

till 

{ai,bi,ci,di}  ♦ {ai,bi,ci,di } . (36) 


The  integral  given  for  the  cross  section  is  then 


N 

,2i  [I- 


hi 


(37) 
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and  our  piece-wise  analytic  expression  for  the  cross  section  is 


d 2a 

dnd|p2l 


a3  |p2|2 

* e2l£il 


N 

,-ri  CI’ 


(38) 


Again  we  stress  that  the  integrations  are  all  performed  analytically.  The  sum 
above  reflects  only  the  subdivision  such  that  W2(x)  and  Wx(x)  can  be 
approximated  with  sufficient  accuracy  by  third  order  polynomials  in  each 


subi nterval . 


Referring  to  the  integrals 


, we  see  that  the  important 


T(l)  t(D  t(2),  t (3)  T (3)  t (4) 
n.i*  !2,i  * Ai  n.i’  2,i*  U.i* 

integrals  to  be  evaluated  are 


and 


m 


dx 


xi-l  Dl1/2(x) 


.p 

!,m  J 


m 


dx 


x i _ i D21/2(x) 


m = -2, -1,0, 1,2  (39) 


= P 

1 ,m  Jv 


m 


dx 


xi-!  Dl3/2(x) 


.p- 


m 


dx 


xi-l  D23/2(x) 


m = -2, -1,0, 1,2  (40) 
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and 


--P  -i 


dx 


X i _1  (2A  - x)  yDj1/2(x)  D?1/2(x), 


(41) 


In  terms  of  these  integrals  we  have 


(i) . 


1 1, i = ai  0oKl,-2  + (bi6o  + aiei)Kl,-l 


+ <Ci60  + bi6l  + ai  62)K1,0 


+ (diS0  + c1e1  ♦ b162)K1>1 


+ C<«i®i  + cl®a>Kl.z  + diB*  3 


(42a) 


1 2 , i " ai  6oK2,-2  + *bi5o  + ai{l*K2,-l 


+ (Cl«0  + + ai«2)K2>0 

+ (di  <S0  + c^!  + b(62)K2jl 


+ (di  5i  + ci 42)K2,2  + di 62  K2,3 


(42b) 


*(i2)  ' a 


1 1 


1 fc'  V b’*nlvTr  Cit><1  " Xi-l)  + Wi(*i  ' X’-l}  (42c) 
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( 42d ) 


I 


I 


I 


I 


I 


(3) 

l,i 


(3) 

2,i 


(4) 

l.i 


(4) 

2,i 


(5) 

1 


aiJl,-2  + biJl,-l 
aiJ2,-2  + biJ2,-l 
aiJl,-l  + biJl,0 

aiJ2,-l  + biJ2,0 


+ ciJ1<0  + d.J1#1 

+ ciJ2.0  + diJ2,l 
+ CiJlj]L  + d-  J^2 
+ ciJ2#1  + diJ2>2 


[ai  + 2Xbi  + (2  A)  2c.  + (2X)  3d.  ] L 

tbi  + 2Xci  + (2X)2di](Jli0  - J2t0) 

[ci  + 2Xdi  ](Jlwl  - J2jl) 


(42e) 

(42f) 

(42g) 


(42h) 


As  just  mentioned,  the  integrals  Ip?  , Ip?  , lj2?  , Ip-  , Ip?  , lj4?  , 
(4)'  (5)' 

\\  , and  Iv  ' follow  from  the  above  expressions  by  substituting  throughout 

, i i 

the  corresponding  parameters  with  primes: 


{^o»^i****}  ' * * * * f 

and 

fai,bi,ci,di}  ■>  (ai,bi,ci,di } . (36) 

Finally,  we  give  the  explicit  expressions  for  the  parameters  which  appear 

in  the  integral  in  terms  of  the  input  variables  for  the  analysis,  viz.,  Ip^l, 

★ 

I p 2 1 , 6,  m,  M,  and  M , as  well  as  the  expressions  for  the  integrals  J,  , 

J.  y 1 1 1 
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Jo  Ki  K0  , and  L in  terms  of  these  variables.  Again  we  leave  the 
2,  m i,m  2,m 

details  of  the  evaluation  to  the  Appendix,  giving  here  only  the  final  results. 
It  is  to  be  noted  that  in  each  case,  the  arqument  of  the  loqarithm  has  been 
written  so  that  there  are  no  cancellations  — all  terms  in  the  argument  are 
positive. 

For  ease  of  writing  we  now  write  and  p2  for  |£1 1 and  |£2I* 
respectively.  Then  with 


e,  = /p  ,2  + m2 


•*2  = 4>22  + m2 


u = 


★ 2 2 

M - M 
2M 


(43) 


(44) 


we  define 


X = £^2  - p1p2cos0  - m2 
|ui|  = (pL2  - 2p1p2cos6  + p22)  1/2 

K = £1  - e2  - W - ft 

2 [X  + (ere2-u)(e1-e2-|_u|)] 

xm  e!-c2- 1 u.1 


2 [X  + (e1.e2-p)(e1-e2+|ju|)] 

XM  " e re 2+| u | 

1 + 


(45) 


(46) 


24 


A,  = 


Bn  = 


eip2’e2picos  9\2  /mp  xsi  n e \2 
p2  + M ) + {— ~) 


- r eip2’e2plcos6' 
2p  ! [ P 1.P  2 - (eiVm  2)cos0]  p2  + 


~2  y 2(  e1P2"  e2P  iCos  0)  + 

4 [p  - y(p1e2-e1p2cos0)  ]2  + 4(ump2sin0)2 


M 


+ -r- 


2m2Pj2(  e1-e2)sin20 


C,  = 


x , = 


>1 


1 A, 


A„  = 


< (e2prelP2cos0)\2  /mp2sin0\2 

vpi  ■ m ) + v- — ; 


Bo  = 2p 2[p  Lp2  - (e^-m^cos©] 


h- 


(47) 


( e2P1” e1P2cos 0)"j  2m2p22(  n20 

+ 


M 


1 


C.  = 


+2yjp1(£2p1-£1p2cos0)  + + 2,qi2)  j 

4 [p2X  + y(p2£1-£2p1cos0)  ]2  + 4(ymp1sin0)2 


M 


x2  " A, 


z,  - 2mp1p2<  sin0 


ni  = f~ 
A1 


^2  = a" 
m2 


*2  v x 7 = «2UX-x2j-  + i)2  1 = m2 


(48) 

(49a) 

(49b) 

(49c) 

2 - 2B2x  + Cj 

(50a) 

2 - 2B2x  + C2 

(50b) 
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B0  = 8m2e1(e1-y)x[p1(p1-p2cos6)  - (e1+e2)u] 


= -2m2 


2Ci(  ery)  ^p1(p1-p2cos0)  - (e1+e2)(icfy)) 


+ X^p1(prp2COS0)  - ( e!+  e2)  vj  ^1  + — j5p) 


2e, 


B2  = m2[p  1(p1-p2cos0)  - (ej+e^ficfyjJfl  + -rp) 


60  = -8m2e2(  e2+p)  X[p2(p2-Plcos0)  + (epe^y] 


<51  = 2m2 


2e2(  e2+w)(p2(P2-Picos0)  + (e1+e2)(,cfw)) 


+ x(p2(p2.plCos0)  + (e1+e2)y)(l  - -~j 


62  = -m 2 [p  2( p 2-p  pos  0)  + ( e!+e2)  ( K+y)  ] (l  - -pp) 
T1  = "X  t*2  " 4m2e1e2  - 2un2(  ef  e2-y)  - 2pe2x] 
t2  = \ lxZ  - 4m2ele2  - 2wn2(  £p£2-y)  + 2y£Lx] 


Pi  = 


JL 

X2 


( X+m 2)  (2 £j e2- X)  + Xk2  + ( X+m2)  ((  X/M)  + y)< 
+ (m2/M2)  X2  + ^X2  + ( X2e1/M) 

( X+m2)  (2  £^2”  X)  + X<2  + ( X+m2)  ((  X/M)  + y)< 


P-,  = 


Y = 


X2  L+  (m2/M2)  X2  + %X2  - (X2e2/M) 

~ [(  X-Hn 2)  (2  e1  e2- X)  + Xk2  + ( X+m2)  ((  X/M)  + y)<] 

A 


n = 1 


(51) 
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30'  = -8m4x[p1(p1-p2cos0)  - (e1+e2)y] 

Sj'  = 4m2[m2(p1(p1-p2cose)  - ( ex+  e2)  ( y)  ) + ^(p  x(  P X-P  2cos  9)  - (e1+e2)y)] 
02'  = -2m2[p1(prp2cos0)  - ( e1+e2)(  ic+y)  ] 

V = 8m4x[p2(p2-Plcos0)  + (e1+e2)y] 

= -4m2[m2(p2(p2-p1cos0)  + ( ( ic+y)  ) + *(p 2(P2"P  icos  e)  + (e1+e2)y)] 

V = 2m2[p2(p2-p1cos0)  + ( e1+e2)(icfy)  ] 


T 


T 


= -2 


fx2  - 2m  2(  X + m2)  1 


2 


X 


. = -21X2  - m4) 

X2 


n'  = -2  . 


(52) 


The  important  integrals  are  then 


°i,o = s9nCAixi  - Bl)  10 


|A1xi  -8,1* 


sgn  fA,x  ■ , - B,  ) £n 

/a7  1 ’-1  1 


^lxi-l  ~ BlJ  * 


(53) 


and 
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, /IB.x,  - c,|  + nhxJ/c.X 

Jl,-1  ■ s«"(Blxi  - Cl)  *>  (— — ) 


/ 1 B, x - , - c.  | + D*(x.  , )/C,  \ 

- “i:  s9n(B1xi_1  - Cj)  In  I " ) (54) 


/C 


xi-le 


where 


sgn(x)  = 


1 

0 

-1 


x > 0 
x = 0 
x < 0 


(55) 


Then  from  q and  ^ we  have 


Jl,l  = t°l(xi)  - Dl(xi-1^  + A^  Jl,0 


(56) 


3B 


1,2  = 2AX  txi  D1  Cxi  ^ + 2Ax  Jl,l  ‘ 2A1  Jl,0 


(57) 


1,-2 


r°iM 


Di  (xi-i } 


1L  Ai  xi-l  J °1 


Si 

+ C7  Jl,-1 


(58) 


Then  the  functions  m are  given  by 


K = 
UO 


(AjX.j-Bj)  ^Alxi-l"Bl  ^ 

°l(xi)  Dlfxi-l) 


(59) 


1,1 


1_ 


(~Bjx  ix  i _ l+C  1 3 

0]_  (x-j  ) ^1  (xi  -1 ) 


(60) 
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(61) 


Kl,2  = Ax  tJl,0  + 2B1K1,1  ’ ClKl,oJ 


K - J- 
Kl,3  Ax 


x - 


i-1 


Dl(xi  _1 ) 


3B, 

+ "A[  Kl,2 


2C 


1 


Aj  Kl,l 


(62) 


Kl,-1  * t 


Dl  (*, ) Dl(xi-l) 


+ C,  Jl,-1  + C,  Kl,0 


(63) 


1,-2 


xiDl(xi) 


3B.  2A, 

+ Cj  Kl,-1  " Kl,0 


(64) 


The  integrals  J0  and  K0  can  all  be  obtained  from  J,  and  K0  merely  by 
2,m  2,m  l,m  2,m 

replacing  all  subscripts  1 by  the  subscript  2 in  the  expressions  above. 

Finally,  we  have  the  integral  L.  In  order  to  write  this  without 
cancellations,  we  define7 


3 = D^(2  X)  = D^(2X)  = 2k/x(X  + 2m2)  (65) 

5l,i  = (2XA1  ‘ Bl)xi  ’ 2XB1  + C1 

h i = <2XA2  " B2^x i " 2XB2  + C2 

?l,i_l  = (2XAx  - B1)xi-1  - 2XBx  + C1 

?2,i_i  = (2xa2  - B2 ) x i-i  - 2xb2  + C2  (66) 


7See  pages  51  and  52  of  the  appendix,  leading  to  eg.  (A. 59). 
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Then 


/|C  | + BD^fx-  )' 


- i s9n(?i,i-i^  ln 


1 ^l.i-1 1 * 


- X sgn(t2j1)  in 


U2  i I + 6D2(xi  ) 


1 r \ 0^2  i-1  ^ + ^xi-l ^ 

+ i sgn(E2>1.1)  *n  (■— ‘-■-1 ‘ L- 


- \ ) - s9n(52,i)jtnlxi  ' 2X1 

+ \ (s9"(el,i-l)  - s9n<52,i-l))*n|xi-l  - 2X1 


(67) 
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APPENDIX 


Algebraic  Details 

In  this  appendix  we  present  some  of  the  detailed  algebra  leading  to  our 
final  expression  for  the  cross  section,  in  particular  to  eqs  (27),  (28),  (38), 
and  ( 44) - ( 67 ) , 

Our  starting  point  is  the  cross  section  integrated  over  the  photon 
azimuthal  angle  given  here  by  eqs  (22)  and  (23),  As  we  have  noted,  this 
is  the  expression  given  by  Mo  and  Tsai  [2],  and  Tsai  [3],  and  Stein  et  al . 

[7],  apart  from  the  minor  modifications  and  typographical  errors  which  we  have 
already  discussed.  Since  we  perform  the  remaining  integration  analytically, 
the  most  appropriate  variable  of  integration  is  q2  rather  than  cose^  in  view 
of  the  functions  W2(-q2)  and  W^-q2).  Further,  as  we  have  already  noted,  the 
form  of  the  integrand  itself  is  then  also  somewhat  simpler.  Referring  to 
eqs  (22),  (23),  and  (24),  we  see  that  our  first  task  is  then  to  express  the 
quantities  listed  in  (24),  as  well  as  the  photon  energy,  u>,  in  terms  of  the 
variable  of  integration,  q2,  and  to  derive  the  differential  change  of  variable 
given  by  (25).  We  see  that  in  (24)  the  variable  cose^  enters  only  either  as  u> 
or  as  cocos 0^.  (The  quantity  b enters  only  as 

b2  = p22sin20pCo2sin20|< 

= p22sin2ep[co2  - (cocosek)2] 

Thus  again  we  have  either  co  or  cocosQ^.) 

Further,  we  note  that  the  expressions  for  cos0p  and  cos0s  listed  in  (24) 
follow  from 
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£i*“  * £i‘(£i  - £2) 


■ P;2  - P]P2C0S6 

■ P]Mcoses  (Pl  = Ipjl)  (A.l) 

and 

£2*U.  - £2*(£j  - £2) 

■ PiP2cos6  - P22 

* P2l£|cosep  (p2*|£2|).  (A. 2) 

In  addition,  it  is  useful  to  note  from  u ■ Pj  - p2  and  the  figure  below, 


that 


! 
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and 


Thus 


p2  px  |u| 


sine  “ sine  " sine  * 
s p 


(A. 4) 


p2sinep  = PjSines 


(A. 5) 


and 


lujsine  = PjSine  (A. 6) 

r 

|u^|sines  = p2sine  . (A. 7) 

The  Jacobian  expressed  by  eq  (25)  follows  directly  from  eqs  (11)  and 
(15) . From  eq  (11),  viz.. 


O) 


★ 2 


u2  - M 
2(u0  - |uJcosek) 


* 


we  have 


(u2  - M*2) |u|d (cose.) 

du>  = = — 

2(u0  - h|cosek)2 

o)|_u|  d(cosek) 

" (u o - |uj  cosek)  * 
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and  from  eq  (15),  viz.. 


w 


el  - 


» 


we  have 


du)  = "2M  d(q2)  * 


so  that 


u)|uj  d (cose.)  , 

(u0  - \u\  cosek)  = 2M  d^q2^ 


or 


2Ma)  d (cos0k)  1 d(q2) 

q4(u0  - |_u|cos0k)  " Tuf  ~q^  * 


which  is  eq  (25). 


Next  we  express  w and  ucos©^  in  terms  of  q2.  In  eqs  ( 12) -( 15)  we  have 
shown  that 


(i)  = £ j 


(A. 8) 


in  which  the  quantities  elte2,u,  and  M are  defined  in  (4).  Next,  from  eq  (6) 
we  have  (four-vectors  are  implied  in  the  next  two  lines) 
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q2  = (pj  - p2  - k)2 


= (P ! - P2)2  - 2k*(p1  - p2) 

= - 2X  - 2^^  - e2)  + 2u>|ujcos0k  (A. 9) 

where,  as  in  (45), 


X = exe2  - p^.£2  - m2  . 


Thus,  using  (A. 8), 

0 - 0 

acos0k  = + (ei-e2Mere2"*0  + %q2Ci  + • (a.io) 

From  (A. 8)  and  (A. 10)  we  note  that  both  co  and  wcos 0^  are  linear  functions  of 
q2.  Thus,  referring  to  (24),  we  see  that  a and  a'  are  also  linear  functions 
of  q2,  and  that  a2  - b2  and  a' 2 - b2  are  quadratic  functions  of  q2. 

The  algebra  involved  in  our  calculation  is  reduced  considerably  by  noting 
that  upon  making  the  substitutions 

subscripts  (1  -*•  2,  2 1),  y - y,  M - M (A. 11) 

we  have,  from  (A. 8)  and  (A. 10), 

to  - 00 
and 

wcos  0^  ♦ cocos©.,  ( i . e . , <cos0k  does  not  change). 
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and  from  (A.l)  and  (A. 2) 


cosQp  ♦ - cos  0^  , cos6s  * " cos0p* 

Thus  when  a and  a'  in  (24)  are  expressed  in  terms  of  q2  as  the  sole  variable, 
we  may  obtain  a'  from  a simply  by  making  the  substitutions  indicated  in  (A. 11) 
and  introducing  an  overall  minus  sign.  That  is,  under  (A. 11) 

a ♦ - a'  , a1  ♦ -a  . 

In  similar  fashion  we  see  from  (A. 5)  that  b2  as  defined  in  (24)  remains 
unchanged  under  the  substitutions  of  (A. 11).  Thus  under  (A. 11)  we  also  have 

a2  - b2  ♦ a'  2 - b2  . 

We  thus  need  only  write  a and  a2  - b2  in  terms  of  q2,  and  then  use  ( A.ll) - 
(A. 13)  to  obtain  a'  and  a' 2 - b2. 

We  start  then  with  a as  given  in  (24),  viz., 

a = e2o)  - p2cos6pwcos9k 

In  this  expression  we  then  substitute  u as  given  by  (A. 8),  cosQp  as  given  by 
(24)  and  wcose^  as  given  by  (A. 10).  After  some  tedious  but  straightforward 
algebra  we  find 


a 


l [p1(p1-p2cos9)  - (e1+e2)M]x 

2 

- + ^q2[p1(p1-p2COS0)  - (e1+e2)(  ic+p)  ] 


(A. 14) 
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Then,  applying  (A. 11)  and  (A. 12)  we  have 


a = - 


[p2(p2-p1cos0)  + (e1+e2)y]X 
+ ^q2[p2(p2-p1cos9)  + (e1+e2)(icfii)] 


(A. 15) 


where  X and  < are  given  by  eg  (43)  in  the  body  of  this  report,  viz., 


X = £^2  - pjPgCosQ  - m: 


K = 


el  - e2  - **  ' M 


(A. 16) 


We  next  consider  a2  - b2,  which,  as  we  have  noted,  is  a quadratic 
function  of  q2.  The  quantity  b is  given  in  (24),  viz., 

b = - t»p  2s  i n 6pS  i n 6^ 

and  from  (A. 6)  this  may  be  written  as 


^lPz  . . 

b = - — |^j|-  sin0sin0k  , 


from  which 


b2  = 


(pxp2sin0) 


[to2 


(ocos0k)2]  . 


We  then  substitute  in  this  expression  u and  ocos0k  as  given  by  (A. 8)  and 
(A. 10),  and  obtain  b2  as  a quadratic  function  of  q2.  With  a as  given  by 
(A. 14)  we  then  have  then  have,  after  much  tedious  but  straightforward  algebra, 
the  following  expression  for  a2  - b2: 
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(A. 17) 


a2  - b2  = {A,  (q 2)  2 + 2B,q2  + C, } 

4u2  1 11 

where  Ap  Bp  and  are  qiven  by  eqs  (47).  It  should  be  noted  from  eqs  (47) 
(in  particular  for  the  integrations  that  follow)  that 

A:  > 0 , C1  > 0 . (A. 18) 


With 


x = -q2  > 0 

as  defined  earlier  in  eq  (26),  we  then  write 

a2  - b2  = — D,(x) 
4u2  1 


(A. 19) 


where 


D^(x)  = A^x2  - 2B^x  + 

= A1[(x  - xx)2  + nx2]  (A. 20) 


with 


x 


l = 


(A. 21) 


and 
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(A. 22) 


V = 


a r . r 2 

Vl  Bl 

A 2 
A1 


With  Ap  B^,  and  as  given  by  (47)  we  find,  after  considerable 

A iC i - Bj2  = (2mp1p2<  sin0)2  = c2  > 0 . 

Again  making  the  substitutions  (A. 11),  we  obtain 


a'  2 - b2  = (A9(q2)  2 + 2B9q2  + C0] 

4u2  2 2 2 


with  A^,  B^,  and  as  given  by  eqs  (48).  We  now  write 


where 


with 


a'  2 - b2  = D9(x) 

4u2  2 


D^(x)  = A^x2  - 2B^x  + 


= A9[(x  - x2)2  + n22] 


x2  = A, 


and 
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algebra, 

(A. 23) 

(A. 24) 
(A. 25) 

(A. 26) 
(A. 27) 


I 


(A. 28) 


From  (A. 23)  we  see  that  A^C^  - B^2 
(A. 11).  Thus 

A2C2  - B22  = A1C1  - Bx2  = c2  . (A. 29) 

Next  we  consider  the  terms  in  (23)  with  factor  v.  From  eqs  (24),  (A. 14), 
and  (A. 15)  we  have 


V?  - B?2 


is  unchanged  under  the  substitutions 


v = 


' I2 

Uq2  + X)  [p2(p2-p1cose)  + p1(p1-p2cose)  ] 
-2 

q2  + 2 X 


(A. 30) 


since,  from  eq  (8b) , 


\j2  = Pi2  _ 2p1p2cos0  + p22  . 

Thus  the  terms  with  factor  v in  the  integral  (22)  have  an  overall  factor 

1 

q 4(q 2 + 2x) 
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As  may  be  seen  by  referring  to  these  terms  in  (23),  they  also  have  a 
coefficient  of  the  form  c0  + c^2  + c2(q2)2.  (For  the  terms  multiplying 
W^-q2)  we  have  cx  = c2  = 0.)  Then,  to  obtain  our  expression  (28),  we  write 


c o + Cjq2  + c2(q2)2 
q“(q2  + 2X) 


c g - 2XCj  + (2X)2 


(2X)2 


(q2  + 2X) 


1 


c0  - 2Xc 
(2X)2 


(A. 31) 


Finally  we  substitute  (A. 14),  (A. 15),  (A. 19),  (A. 25),  and  (A. 30)  with 
x = - q2,  in  the  cross  section  given  by  (22)  and  (23),  using  (A. 31)  for  the 
terms  with  factor  v.  Collecting  the  terms  with  similar  denominators  then 
gives  the  expression  (28). 

We  derive  next  the  expression  given  in  (46)  for  the  limits  of 
integration,  xm  and  x^.  From  (11),  viz.. 


we  see  that  <*>  is  an  increasing  function  of  cose^,  and  from  (15),  viz.. 


it  follows  that  q2  is  also  an  increasing  function  of  cos0^,  and  hence  that 
x = -q2  is  a decreasing  function  of  cose^.  Thus  the  minimum  value  of  x, 
viz. , 


o> 


(ID 


q2 

“ = el  - c2  - M + 2M  • 


(15) 


corresponds  to  cos0k  = + 1 
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and  the  maximum  value  of  x,  viz.. 


x = x^  = -qM2  corresponds  to  cose^  = - 1 


We  now  compare  eq  (15)  given  just  above  with  eq  (A. 10),  viz., 


G — G 

“COS  9.  = -At[x  + (e1-e2)(£1-c2-u)  + *q2(l  + 


When  cos 9^  = +1  these  two  equations  give 


’m 


2 1 


el-e2-y  + 2M  " | u | 


[x  + ( Ej- e2)  ( ei"  e2"y)  + + m 


e , - 1 


or 


2[X  + (e1-e2-|u>|)(ere2-|i)] 


m 


1 + 


WW 


as  in  (46).  When  cos0k  = -1  these  two  equations  give 


el’  e2-li  + 2M 


?m 7T7T[x  + + Hi^1 


or 


XM  ' 


2 [X  + ( er£2+|_u|)(  ere2-u)  ] 


1 + 


er  e2+  'J£' 


as  in  (46). 


(A. 10) 


-)] 
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We  next  show  that  the  point  x = -q2  = 2A,  where,  as  in  (45), 

A = £^2  - p xp 2cos  0 - m2  , 

lies  within  the  range  of  integration: 

x < 2A  < xM  , 
m M 

as  noted  in  (29).  Now  from  (A. 9)  we  see  that  x = -q2  = 2 A implies 

el"  e2 

cos0k  = |ui|  * (A. 32) 

Now  as  we  show  immediately  below, 

£ - G 

0 < -ijy2  < 1 . (A. 33) 

Thus  since  the  integration  region,  when  we  choose  cos0k  as  integration 
variable,  is  the  entire  range  of  possible  values,  - 1 < cos0k  < 1,  we  see  from 
(A. 32)  and  (A. 33)  that  the  point  x = 2A  corresponds  to  a point  in  the  range  of 
integration  with  x as  integration  variable  (recall  that  x is  a uniformly 
decreasing  function  of  cos0k),  i.e.,  that 

x < 2A  < xM  . 
m M 

Further,  since  the  photon  energy,  given  by  (15),  must  be  positive  for  cosek 
within  the  range  of  integration,  i.e.. 
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a)  = e,  - e 


- y + > 0 , 


1 2 - - ■ 2M 


it  follows  that  this  is  also  true  for  the  point  - q2  = 2X  just  considered. 
Thus 


<=ei"e2"w"M>0  * 


(A. 34) 


Returning  now  to  the  proof  of  (A. 33),  we  note  from  e2  < e1  (and  p2  < p^  that 


:l  ’ pl  = e,  + 


m2  „ m2 


l + Pi  e2  + P2 


= e„  - 


2 - P2 


Thus 


0 < ere2  < prp2  < |prp2l 


from  which  (A. 33)  follows. 

Finally,  we  conclude  our  presentation  of  algebraic  details  with  a 
consideration  of  the  expressions  given  in  eqs  (53) - ( 67)  for  the  integrals 
defined  in  eqs  ( 39) - ( 41 ) . We  begin  with  g,  defined  by  eq  (39),  viz.. 


wf 


dx 


xM  D1l/2(x) 


(A. 35) 


where  from  eq  (50a) 


D^(x)  = A^x2  - 2BjX  + C ^ 
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A straighforward  integration  gives  [9] 


_L  £n  /Alxi-1  " B1  * D1  (xi-l)v^l 


/A. 


(A. 36) 


where  c0  is  an  arbitrary  constant.  Now  for  reasons  that  will  be  apparent 
momentarily,  we  choose 


c o = AjCj  - B^2 

= ? = 2mp1p2ic  sin0  > 0 (A. 37) 

from  (A. 23)  and  (A. 34).  The  arguments  of  the  logarithms  in  (A. 36)  are  then 
dimensionless  and,  from  eqs.  (47)  and  (A. 23),  are  of  order  e/m  if  A^x  - > 0 

(where  x = x.  or  x-_j).  0n  the  other  hand,  since 

D1(x)A1  - (Axx  - Bx)2  = (Axx2  - 2Bxx  + C^Aj  - (AjX  - B^2 

= A1C1  - Bx2  (A. 38) 


we  have 


Axx  - B1  + D1^(x)/A1 


-(Ajx  - Bx)  + D1?5(x)/A1 


(A. 39) 
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Thus  if  A^x  - < 0 we  can  write  the  arguments  of  the  logarithms  in  (A. 36)  in 

the  form  given  by  the  right  hand  side  of  (A. 39),  from  which  we  see  that  the 
argument  is  then  of  order  m/e.  Thus  if  A^x  - < 0 there  is  a cancellation, 

leaving  terms  of  relative  order  (m/e)2.  We  will  therefore  write  the  argument 
so  that  there  is  no  cancellation,  i.e.,  as  given  by  the  left  hand  side  of 
(A. 39)  if  AjX  - Bj  > 0 and  by  the  right  hand  side  of  (A. 39)  if  A^x  - B^  < 0. 

These  can  be  combined  into  a single  expression  by  writing 


|Alxi  - Bi!  * 


1 

— sgn(A1xi.1  - Bj)  An 


1 ' i-l;  1 


/A1 


where 


1 


x > 0 


sgn(x)  = 0 x = 0 . 


-1 


x < 0 


This  is  the  result  given  in  (53). 

The  integral  ^ is  also  defined  by  (39),  viz.. 


(A. 40) 


Again,  a straightforward  integration  gives  [9] 
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+ — ln  rBixi-i  + ci  + Di^xj-i)/ci 


/c 


(A. 41) 


1 


xi-lco' 


where  c0'  is  an  arbitrary  constant.  In  analogy  with  (A. 38)  we  note  that 


D1(x)C1  - (-BlX  + C^2  = (AjX2  - 2Bxx  + - (-B^  + C^2 


= x2(A1C1  - Bx2) 


= x2?2  . 


(A. 42) 


Thus  we  again  choose  c0'  = c,  and  note  that 


-B^x  + C-^  + D-j^(x) 


xc 


x C 


-(-BjX  - Cj)  + DjVj/Cj 


(A. 43) 


Thus  if  -B^x  + Cj  < 0 we  can  write  the  arguments  of  the  logarithms  in  (A. 41) 
in  the  form  given  by  the  right  hand  side  of  (A. 43),  thus  avoiding  the  cancel- 
lation. As  before,  we  now  write  ^ in  a form  which  avoids  the  possible 


cancellation  in  the  arguments  of  the  logarithms,  viz.. 


1 /l-B.x,  + C,|  + D.^xJ/C, 

Ji-i ' ' *=  ‘•"'•V*  * c.’  ” ( 1 , . 


1 , + C.  | + D,  ^(x . ,)/C, 

+ ^ *9nC-81*1.1  ♦ Cx)  - ^ 1 ^ 


In 


xm' 


which  is  equivalent  to  the  result  given  in  (54). 

The  integrals  g and  defined  by  (40)  and  given  expl icitely  by 

(59)  and  (60),  are  elementary  and  present  no  particular  difficulty.  The 
expression  of  m and  m in  terms  of  g,  o*  and  K1  1*  9^ven 

( 56 ) - ( 64) , is  also  elementary  and  may  be  found,  for  example,  in  the  integral 
tables  of  Grffbner  and  Hofreiter  [9]. 

Lastly,  we  have  the  integral  L,  defined  in  (41),  viz., 


L = 


1 

(2X  - x) 


.D1l/2(x)  D21/2(x) 


dx 


(A. 44) 


We  first  note  that 


D1(2  X)  = D2(2  X)  , (A. 45) 

and  thus  the  integrand  has  in  fact  no  singularity  at  x = 2X.  This  may  be  seen 
most  easily  as  follows:  From  (A. 30) 

a'  - a = -^(q2  + 2X) 


= Mx  - 2 X) 


(A. 46) 
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and  from  (A. 19)  and  (A. 25) 


Thus 


and 


We  define8 


Dx(x)  - D2(x)  = 4u_2 (a 2 - a'2)  . 


Dx(x)  - D2(x)  = -2jj2(a  + a')(x  - 2\) 


Dx(2X)  = D2(2X) 


3 = D11/2(2X)  = D21/2(2X)  . 

In  analogy  with  the  procedure  followed  for  q and  ^ we  first 
integration  straightforwardly  (see  [9])  and  write 


L = 


i In  N.i  * BDltxi] 

6 *"  \ |xr2X|c 


1 tn  /sl,j-l  + gDitxj-i) 
» \ |xi-l_2X|  5 , 


1 >n 

S \ |xi-2X|c  , 


8See  the  end  of  this  appendix  for  a derivation  of  the  value  of  3. 


(A. 47) 

(A. 48) 

(A. 49) 
perform  the 


(A. 50) 
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where 


^(x)  = (2 XA1  - B t )x  - 2XBx  + Cl 


C2(x)  = (2 AA2  - B2)x  - 2XB2  + C2 


(A. 51) 


5l,i  ~ Mxi^  * 5l,i-l  = Mxi-1^ 


^2,i  = Vxi^ 


^2,i-l  = ^2  ^ x i - 1 ^ 


(A. 52) 


Then  noting  that 


^(x)  + BD^(x  ) I x -2 X | ? 

|x-2A|c  -^(x)  + SO^(x) 


(A. 53) 


(and  a similar  equation  in  which  all  subscripts  1 are  replaced  by  the 
subscript  2),  we  start  by  writing  L in  a form  similar  to  that  just  used  for 
q and  so  that  there  is  no  cancellation  when  any  of  the  terms  ^ 

52,i-  or  52,i-l  is  native: 


1 ( \ 1 1 ^l.i  I + (xj  ) 

L = 6 sgn(^l,i)  ln  \ | x 1 -2 X | c 

1 /ISi  i J + eofy.  ,) 

- 1 *n  ( ' ix^^-zxu  ' 

" B sgn(^2,i  ) *n  \ | x - -2  X | c j 


i /1 , 1 1 + eoJ(x.  , ) 

+ 1 ^2,1.1) *n  ( tT-pFirr^- 


(A. 54) 
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However,  in  this  form  of  the  expression  for  L there  can  be  a large 
cancellation  for  or  x-  ^ sufficiently  close  to  2A,  since,  as  we  have  noted, 
the  integrand  in  the  integral  for  L,  (A. 44),  has  no  singularity  for  x = 2A. 

We  therefore  write  L in  a form  such  that  this  cancellation  is  manifest. 
Equation  (A. 54)  may  simply  be  written  in  the  form 


and  a similar  expression  with  i replaced  by  i-1.  This  is  the  result  given  in 
eq  (67).  Now  from  (A. 51)  and  eqs  (50a)  and  (50b)  we  have 


(A. 55) 


where 


(A. 56) 


5:(2X)  = Dx(2X) 


S2(2\)  = D2(2A) 


(A. 57) 


and  hence  from  (A. 48)  and  eqs  (50a)  and  (50b) 
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^(2X)  = ^(2\)  > 0 


(A. 58) 


Therefore,  for  x.  sufficiently  close  to  2X  we  have  sgn( 5.  .)  = +1, 

t * 


sqn(£?  ,)  * +1  and  thus  L.  in  (A. 56)  can  be  written  more  simply  as 


Here  (i.e.,  when  sgn(  5-,  J = sgn(C7  • ))  there  is  clearly  no  singularity  for 
x-  = 2X.  These  observations  can  be  used  to  avoid  any  possible  problem  arising 
in  a computer  program  from  x^  or  x-_^  being  chosen  too  close  to  2X.  We  may 
simply  define  L by  (A. 55),  with  L-  (and  L-^)  given  by  (A. 59)  if  both  ^ . > 0 
and  . > 0,  and  otherwise  given  by  (A. 56).  (It  can  be  shown  that  we  can  not 


have  both  ■ < 0 and  C?  1 < 0.) 

Finally,  we  require  the  value  of  6,  defined  in  (A. 49),  given  in  the  body 
of  the  report  in  eq  (65),  and  appearing  in  the  equations  for  L,  (A. 50)  and 
(A.54)-(A.59) . From  the  definition,  (A. 49), 


it  is  clear  that  one  can  evaluate  6 simply  by  substituting  x = 2X  in  the 
equations  defining  D^(x)  and  D^x),  (50a),  and  (50b).  This  is  a rather 
cumbersome  procedure,  however.  We  therefore  give  the  following  simpler 
derivation: 

For  x = 2X  we  have 


(A. 59) 


6 = dJ(2A)  = d|{2X)  , 


a 


a 


D1(2X)  = D2(2X)  , 
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from  (A. 46)  and  (A. 47),  and 


(D  = < 


from  eqs  (15)  and  (A. 16).  Further,  also  for  x = 2X,  we  have 


cos0k  = 


from  (A. 32),  from  which 


sin6k 


/u2  - '(ere2)a 

lul 


✓2X 
' 111 


(A. 60) 


Substituting  these  values  for  w,  cose^,  and  sine^  in  the  expressions  (24)  for 
a'  and  b we  have,  for  x = 2X, 


a'  = ^ Ieiii2  - Pi(e1-e2)(PrP2cos9)  ] (A. 61) 

u2 


and 


<p  ,p2  si  n0  /2X 

b = - . (A. 62) 

u2 

In  writing  (A. 61)  we  have  used  the  expression  for  cos6s  given  in  (24),  and  in 
writing  (A. 62)  we  have  used  the  expression  for  sin0p  given  in  (A. 6). 

Recalling  that  u2  = px2  - 2p1p2cos6  + p22,  a little  algebra  then  gives 
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a'  = — (e1+e2)(e1e2  - p^cosQ  - m2) 


<X(ei+e2) 


(A. 63) 


Substituting  (A. 63)  and  (A. 62)  in  (A. 25)  we  then  have 


D2 (2 X)  = [X(e1+e2)2  - 2p  12p22si n2e]  . 

^ 1 1 2 


Again  recalling  that  X = - p^cosQ  - m2,  a little  algebra  gives 


Mei+e2)2  * 2p12p22sin20  = (exe2  - Plp2  + m2)^2  - 2p1p2cos6  + p22) 


so  that 


= ( X + 2m2)  u2  , 


D 2 ( 2 X)  = 4 k2X(  X + 2m2) 


(A. 64) 


and 


e = D^(2 X)  = D^(2X)  = 2</x(  X + 2m2)  . 


(A. 65) 
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The  form  factors  used  in  these  timing  runs  were  derived  from  a DWBA 
calculation  for  12C  at  500  MeV  (see  subroutine  PWAFFI  in  the  computer 
program  listing,  pp.  98-100). 

PWA  denotes  the  piece-wise  analytic  evaluation  of  the  integral  as 
presented  in  this  report. 

All  times  listed  are  in  seconds.  Numerical  integrals  were 
calculated  with  an  accuracy  of  0.0001 H. 

p ^250  MeV 


1 deg . 


P2  (MeV) 

Simpson 

Gauss  <n=16) 

Kronrod  (n:61) 

PWA 

10 

21.07 

5.42 

1.57 

. 43 

125 

62.03 

16. 12 

2.36 

.22 

240 

3.15 

1.22 

. 05 

.02 

5 deg . 

P2  (MeV) 

Simpson 

Gauss  (n:16) 

Kronrod  <n=61) 

PWA 

10 

11.66 

2.06 

.96 

.42 

125 

4.37 

1.30 

.96 

.22 

240 

1.86 

.59 

.61 

.03 

45  deg. 

P2  (MeV) 

Simpson 

Gauss  <n=16) 

Kronrod  (n=61) 

PWA 

10 

2.46 

.84 

. 87 

.43 

125 

2.33 

1.00 

.78 

.23 

240 

2.21 

.07 

.80 

.02 

90  deg. 

P2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  (n=61) 

PWA 

10 

3.01 

1.24 

. 93 

.44 

125 

2.61 

1.17 

.84 

.23 

240 

2.46 

1.23 

.09 

.02 

135  deg. 

P2  (MeV) 

Simpson 

Gauss  (n:16) 

Kronrod  <n=61) 

PWA 

10 

5.35 

1.78 

1.17 

.43 

125 

3.15 

1.37 

.92 

.23 

240 

3.30 

1.31 

.97 

.01 

175  deg. 

P2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  (n=61) 

PWA 

10 

11.71 

2.43 

1.06 

. 44 

125 

4.47 

1.73 

1 . 18 

.23 
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All  times  are  in  seconds.  Numerical  integrals  uiere 
calculated  with  an  accuracy  of  0.0001%. 


Pl=500  MeV 


1 deg . 

p2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  (n:61) 

PWA 

10 

48.45 

8.53 

2.39 

.99 

250 

70.24 

16.29 

2.90 

. 44 

490 

3.66 

1.20 

1.01 

.03 

5 deg. 

P2  (MeV) 

Simpson 

Gauss  <n=16) 

Kronrod  (n=61) 

PWA 

10 

33. 19 

3.06 

2.21 

.99 

250 

6.57 

2.16 

1.32 

. 44 

490 

3.35 

1.24 

1.04 

.03 

45  deg. 

P2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  <n=61) 

PWA 

10 

3.26 

1.31 

.75 

.98 

250 

3.76 

1.64 

1.33 

.46 

490 

5.01 

2.32 

1.96 

. 02 

90  deg. 


P2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  (n=61) 

PWA 

10 

6.40 

2.33 

1.29 

.99 

250 

4.31 

1.75 

1.55 

.48 

135  deg. 

P2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  (n=61) 

PWA 

10 

20.92 

3.61 

1.97 

.99 

250 

6.40 

2.58 

1.44 

.53 

175  deg. 

P2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  (n=61) 

PWA 

10 

23.49 

5.90 

1.86 

1.00 

250 

8.70 

3.46 

1.89 

.56 
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All  times  listed  ere  in  seconds.  Numerical  integrals  were 
calculated  with  an  accuracy  of  0.0001*. 

Pj-750  MeV 


1 deg . 


P2  (MeV) 

Simpson 

Gauss  <n=16) 

Kronrod  <n=61) 

PWA 

10 

93.26 

16.85 

2.52 

1.63 

375 

46.09 

18.82 

2.96 

.70 

710 

5.72 

1.86 

1.23 

.08 

5 deg. 

P2  (MeV) 

Simpson 

Gauss  ( n : 16 ) 

Kronrod  <n:61) 

PWA 

10 

59.96 

15.38 

3.94 

1.63 

375 

7.26 

2.33 

1.36 

.70 

710 

4.49 

1.89 

1.27 

.08 

45  deg. 

P2  (MeV)  . 

Simpson 

Gauss  (n=16) 

Kronrod  <ns61) 

PWA 

10 

5.66 

2.33 

1.02 

1.62 

375 

5.00 

2.09 

1.75 

.75 

710 

5.60 

2.57 

1.30 

.07 

90  deg. 

P2  (MeV) 

Simpson 

Gauss  (n=16) 

Kronrod  (n=61) 

PWA 

10 

11.47 

4.40 

1.92 

1.62 

375 

6.68 

3.44 

1.94 

1.21 

135  deg. 

P2  (MeV) 

Simpson 

Gauss  (n:16) 

Kronrod  (n=61) 

PWA 

10 

26. 18 

7.16 

2.36 

1.62 

375 

9.84 

2.39 

2.14 

.98 

175  deg. 

P2  (MeV) 

S i mpson 

Gauss  (n=16) 

Kronrod  <n=61) 

PWA 

10 

57.  47 

15. 11 

3.70 

1.62 

375 

10.22 

3.56 

1.47 

.88 
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COMPUTER  PROGRAM  LISTING 


In  the  pages  which  follow  we  give  the  listing  of  the  computer  program. 
However,  the  subroutines  PROMPT,  TOD,  and  CLK  (which  are  called  in  the  main 
program  MSPLTIM,  as  well  as  in  the  timing  subroutines  STARTT  and  STOPT)  are 
system  information  routines  which  are  unique  to  the  Perkin-Elmer  system  on 
which  the  program  was  run.  They  are  not,  therefore,  included  in  the  program 
listing. 


I 


1 

STITL 

MSPLTIM  — timing  for  Tsai  and  PWA 

2 

' C 

This  program  compares  the  formula  for  the  radiative  tail 

3 

C 

given  by  Tsai  and  the  Piece-Wise  Analytic  method 

( PWA ) 

4 

C 

5 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

6 

c 

Maximum  number  of  angles 

• 

7 

PARAMETER  (MNTslO) 

e 

c 

Maximum  number  of  momenta 

9 

PARAMETER  (MNP=10) 

10 

LOGICAL  OPENED, FAST 

n 

DOUBLE  PRECISION  TIM( 4 ) , THETA( MNT ) , PF IN ( MNP ) 

12 

INTEGER  ICLKBC4) 

13 

CHARACTER  FILE*40, T IME*8 , DATE*8 , YESNO*SO 

14 

COMMON  /EUICOM/  NSTEP,EPS 

15 

COMMON  /'CONST / EM,  EM2,  AM,  HBARC,  ALPH,  PIE,  ALPHA,  RTD 

, DTR 

16 

COMMON  /KIN/  EEX, TH, AMO, S3, P3, CTH, STH, S32, ES, P32, 

EP, AMF,SDP3, SDP, 

17 

+ U32, U3, UO, U2, CTHS, CTHP, STHP, XMU 

18 

COMMON  /CORREC/IWTILD, WTIL1, WTIL2, IS0FT,S0FT1 

19 

COMMON  /TIME/  T I M, I CLKB 

20 

DATA  FAST/. FALSE./ 

21 

CALL  TOD ( T IME , DATE ) 

22 

CALL  INICON 

23 

INQUIRE  <1,0PENED=0PENED) 

24 

IF  ( .NOT. OPENED)  OPEN  ( 1 , FILEs  * CON : * ) 

25 

2 

IF  ( .NOT. OPENED)  THEN 

26 

WRITE  (1,901)  ' ', TIME, DATE 

27 

901 

FORMAT  (A, 'MSPLTIM  version  1A  ',A,1X,A, 

28 

+ ' (last  edit  09/06/83)') 

29 

CALL  PROMPT( 1, 'Output  file') 

30 

END  I F 

31 

READ  (1,80)  FILE 

32 

IF  ( FILE . EQ . ' ')  STOP. 'All  done' 

33 

00 

FORMAT  (A) 

34 

OPEN  (6, FILErFILE, ACCESS =' SEQUENT  I AL ' ) 

35 

IF  ( .NOT. OPENED)  CALL  PROMPT ( 1 ,' Exc i t at i on  energy 

' ) 

36 

READ  (1,*)  EEX 

37 

EEXrEEX/HBARC 

38 

IF  ( .NOT. OPENED)  CALL  PROMPT ( 1 ,' Inc i dent  momentum 

' ) 

39 

READ  (1,*)  S3 

40 

S3: S3/HBARC 

41 

IF  ( .NOT. OPENED)  CALL  PR0MPT(1,'A  of  target') 

42 

READ  (1,*)  AMO 

43 

„• 

AMO: AMO*AM 

44 

IF  ( .NOT. OPENED)  WRITE  (1,910) 

45 

910 

FORMAT  ('  Enter  values  for  final  momentum  (end  with  negative):') 

46 

NP  : 0 

47 

99 

READ  (1,*)  P3 

48 

IF  (P3. GE. 0. ) THEN 

49 

NP : NP+ 1 

50 

PFIN(NP) :P3/HBARC 

51 

IF  (NP.EQ.MNP)  GO  TO  88 

52 

GO  TO  99 

53 

END  I F 

54 

IF  (NP.EQ.O)  STOP 

55 

88 

IF  ( .NOT. OPENED)  WRITE  (1,909) 

56 

909 

FORMAT  ('  Enter  values  for  THETA  (end  with  negative):') 

57 

NT  : 0 

60 


5B 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

103 

109 

110 

1 1 1 

112 

113 

,114 


.9  READ  (1,*)  TH 

IF  (TH. GE. 0. ) THEN 
NTrNT+1 

THETA ( NT )=TH*DTR 
IF  (NT.EQ.HNT)  GO  TO  9 
GO  TO  8 
ENDIF 

IF  (NT.EO.O)  STOP 

9 IF  ( .NOT. OPENED)  CALL  PROMPT ( 1 , 

+ 'Starting  number  of  steps  for  the  integral') 

READ  (1,*)  NSTEP 

IF  ( .NOT. OPENED)  CALL  PROMPTt 1, 'Accur acy  of  the  integral') 
READ  (1,*)  ACC 

IF  ( .NOT. OPENED)  CALL  PROMPT ( 1 , 'W-TILDA  correction') 

READ  (1/80)  YESNO 
IF  ( INDX ( YESNO/ ' Y ' ) .GT.O)  THEN 
IWTILD=1 

ELSE 

IWTILDrO 
ENDIF 
ISOFTrO 
CALL  PWAFFI 
DO  502  1=1, NT 

TH=  THETA ( I ) 

WRITE  (6/901)  ' 1 TIME/ DATE 
WRITE  (6/902)  TH*RTD, S3*HBARC , 

+ AMO/AM/ EEX*HBARC/ ACC, YESNO 

902  FORMAT  ('0  Scattering  angle  =',F7.2,'  degrees'/ 

+ T62/  ' Inc  ident  momentum  =',F8.3, ' MeO/c'/ 

+ ' Mass  of  nucleus  ='/F8.3/'  AMU'/ 

+ T61 , ' Exc i t at i on  energy  = ',F8.3,'  MeU'/ 

+ ' Accuracy  of  numeric  integral  =',1PE8.1, 

+ T6l/'W-Ti)da  correction  made?  '/A5) 

WRITE  (6/903) 

903  FORMAT  ( ' 0 ' , T6 . 'F i n a 1 ' / T17,  ' Ts a i ' / T26/ ' For m factor'/ 

+ T39, ' Integr al ' , T52/ 'Total '/T65, 'PWA',T76/ 

+ 'Search '/ T37/ 'Integral '»T100/ 'Total ',T110/ '*  Diff' 

+ T126/ 'Speed ' / 

+ T3/'momentum'/T15/'result'/T27/ 'time'/ T39/ 'time'/ 

+ T51 / 'time', T63, ' resu 1 t ' , T75, 'time', T87, 'time', 

+ T99, ' t i me ' , T124,  ' incre ase  ' ) 

CALL  CLK(ICLK,Y'2B', .TRUE. ) 

DO  503  J=1,NP 
P3=PF IN ( J ) 

C Initialize  some  commons 

CALL  INIKIN 
T I M ( 1 )=0. DO 
TIM(2)=0. DO 
TIM(3)=0. DO 
TIM(4)=0.D0 

IF  ( J . LE . 1 . OR . . NOT . FAST ) THEN 
EPS=-ACC 

ELSE 

EPS=ACC*SIGTS 

ENDIF 

CALL  ST  ARTT ( 2 ) 

CALL  TSAINT (SIGTS) 
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115 

CALL  STOPT ( 2 ) 

116 

• 

CALL  ST ARTT ( 4 ) 

117 

CALL  PWASPL(SIGPWA) 

118 

CALL  STOPT ( 4 ) 

119 

WRITE  (6/905)  P3*HBARC  *SIGTS»TIM(1)»TIM(2>-TIM(1)» 

120 

♦ 

TIM(2)»SIGPWA»TIM(3)*TIM(4)-TIM(3)»TIM(4)> 

121 

+ 

200.  *<SIGPWA-SIGTS>/'<SIGPWA  + SIGTS) , T I M < 2 ) ✓T IM  ( 4 ) 

122 

905 

FORMAT  (1X#F8.2*T12» 1 PE  1 2 . 5 » T24 , 0PF9 . 3 , T36, F9. 3,T4e, 

123 

♦ 

F9. 3,T60, 1PE12. 5/T72, 0PF9.3, T84,F9. 3, T96/F9.3, T108, 

124 

+ 

1PE10. 3/ T 1 20 » 0PF10. 2) 

125 

503 

CONTINUE 

126 

502 

CONTINUE 

127 

CLOSE 

(6) 

128 

GO  TO 

2 

129 

END 
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6 

7 

8 
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• $TI TL  STARTT  --  Start  a tune  interval 
SUBROUTINE  STARTT(I) 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 
DOUBLE  PRECISION  T I M f 4 ) 

INTEGER  I CLKB ( 4 ) 

COMMON . /TIME/  TIM, ICLKB 

CALL  CLK(ICLKB(I),Y'2B', .FALSE. ) 

RETURN 

END 
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7 
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ST1TL  STOPT  — Stop  a time  interval 
SUBROUTINE  STOPT(I) 

C Note  that  the  constant  is  for  optimized  Fortran. 
IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

DOUBLE  PRECISION  TIM<4) 

INTEGER  ICLKB ( 4 ) 

COMMON  /TIME/  TIM, ICLKB 
CALL  CLKdCLKE,  Y'20',  .FALSE.  ) 

TIMC I)=TIM(I )+( ICLKE-ICLKBC I ))/ 1000 0009374 

RETURN 

END 
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STITL  INICON  — Initialize  constants 
SUBROUTINE  INICON 

C Subroutine  to  initialize  the  common  CONST.  All  energies  in 

C units  of  fm-l.  Constants  from  RMP  52*  No.  2 (April  1980)  as 

C reprinted  by  the  CERN  particle  properties  data  booklet. 

IMPLICIT  DOUBLE  PRECISION  (A-H*0-Z> 

COMMON  /CONST/  EM* EM2* AM* HBARC* ALPH* PIE* ALPHA, RTD* DTR 

HBARC: 197. 32858D0 

ALPH=137. 03604D0 

AMr 931 . 5016D0/HBARC 

EM= . 51 10034D0/HBARC 

EM2:EM**2 

PIE=3. 1415926535897D0 
ALPHA= 1 .ODO/ALPH 
DTR= PIE/ 180. DO 
RTD:180. DO/PIE 
RETURN 
END 


65 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 


no  on 


$TITL  INIKIN  — Initialize  kinematic  variables 
SUBROUTINE  INIKIN 

C Subroutine  to  initialize  the  commons  KIN  and  CORREC. 

C In  general  this  must  be  called  if  any  of  the  first 

C five  elements  of  KIN  change. 

IMPLICIT  DOUBLE  PRECISION  (A-H/O-Z) 

COMMON  /CONST/'  EM,  EM2,  AM,  HBARC,  ALPH,  PIE/ ALPHA,  RTD, DTR 

COMMON  /KIN/  EEX, TH, AMO, S3, P3, CTH, STH, S32, ES, P32 , EP , AMF , SDP3 , SDP. 

+ U32,U3,U0,U2,CTHS,CTHP,  STHP,XMU 
COMMON  /CORREC/IWTILD, WTIL1, WTIL2, IS0FT,S0FT1 

C sin's  and  cosin's  of  angles 

CTHrDCOS(TH) 

STH=DSIN(TH) 

C Incident  momentum  squared 

S32=S3*S3 

C Incident  energy  of  the  electron 

ES=DSQRT ( S32+EM2 ) 

C Final  momentum  squared 

P32=P3*P3 

C Final  energy  of  electron 

EPrDSQRT ( P32+EM2 ) 

C Final  rest  mass  of  the  nucleus 

AMF  = AMO+EEX 

C This  is  ( AMF**2  - AM0**2 ) /AMO 

XMU=EEX+. 5D0*EEX**2/AM0 

C 3-vector  dot  product  of  s and  p 

SDP3=S3*P3*CTH 

C 4-vector  dot  product  of  s and  p 

SDP  =ES*EP-SDP3 

C Momentum  part  of  u=s+t-p  (momentum  loss  of  the  electron) 

U32=S32+P32-2.D0*SDP3 
U3=DSQRT(U32) 

C Energy  part  of  u=s+t-p 

UO=ES+AMO-EP 

C Length  of  u=s+t-p  (squared)  (missing  mass  squared) 

U2=U0**2-U32 

C cosin  of  angle  between  s and  u 

CTHS=(S3-P3*CTH)/U3 

C cosin  of  angle  between  p and  u 

CTHP=(S3*CTH-P3 ) /U3 

C sin  of  angle  between  p and  u 

STHP=S3*STH/U3 

Initialize  coefficients  for  change  of  W to  W-TILDA. 

WTIL1=1 . D0+28.D0*ALPHA/(9.D0*PIE)-ALPHA*DL0G(ES/EP)**2/(2. DO*PIE)+ 
+ ALPHA* (PIE*P IE/6. DO-SPENCE( . 5D0*( 1 . DO+CTH) ) )/PIE 
WTIL2=13. D0*ALPHA/(PIE*6. DO) 

Initialize  coefficients  for  soft  photon  correction. 

WS  = ES-EP/( 1 . DO— EP* ( 1 . DO-C  TH ) /AMO ) 

WPr ES/ ( 1 . DO+ES* ( 1 . DO-CTH ) /AMO ) -EP 
SOFTlrWS*WP/(ES*(EP+WP) ) 

RETURN 
END 
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STITL  TSAINT — Evaluate  Tsai  by  doing  numerical  integral 
SUBROUTINE  TSAINT ( UAL ) 

C The  formula  given  by  Tsai  in  Slac-pub  848  (as  we  1 1 as  in 

C Mo  and  Tsai  RMP  41,  205  and  Stein  et . al.  Phys  Rev  D 12,1884) 

C is  evaluated  by  doing  the  integral  numerically.  This 

C is  done  in  up  to  7 separate  intervals.  Cal.l  INIKIN  before 

C calling  this  routine. 

IMPLICIT  DOUBLE  PRECISION  CA-H,0-Z> 

EXTERNAL  TSAI 
LOGICAL  DEBUG 
INTEGER  NS ( 7 ) 

DOUBLE  PRECISION  EDGE ( 8 ) 

COMMON  /CONST/  EM,  EM2,  AM,  HBARC, ALPH, PIE, ALPHA, RTD, DTR 
COMMON  /KIN/  EEX , TH, AMO, S3, P3, CTH, STH, S32,  ES, P32, EP, AMF, SDP3, SDP, 
+ U32,U3,U0,U2,CTHS, CTHP, STHP,XMU 
DATA  DEBUG/. FALSE. / 

DATA  EDGE C 1 ) /-I . DO/, EDGE C 8 ) /I . DO/ 

UAL=0. DO 

THL IM=AM0*CES-EP)/C2. DO*ES*EP) 

IF  (THLIM.LE. 1. DO)  THEN 

THLIM=2. DO*DASIN( DSQRT(THLIM) ) 

IF  (TH.GT.THLIM)  RETURN 
ENDIF 

THS  = DATAN2  C DSQRT  C 1 . D0-CTHS**2) , CTHS) 

DTHS= DSQRT CEM/ES)*2.0D0 
THP=DATAN2 (STHP, CTHP) 

DTHP=DSQRTCEM/EP)*2.0D0 

EDGE (2) =DC0S(DMIN1 (PIE, THP+DTHP) ) 

EDGE ( 3 ) =CTHP 

EDGE (4) rDCOSCDMAXl (0. DO, THP-DTHP ) ) 

EDGE (5) = DCOS  C DMIN1 CPIE, THS  + DTHS)  > 

EDGE ( 6 ) = CTHS 

EDGE (7 ) = DCOS  C DMAX1 (0. DO, THS-DTHS) ) 

IF  CEDGEC4).GT.EDGEC5) ) THEN 

EDGE C 4 ) = CEDGEC3)+EDGE(6) )/2. DO 
EDGEC5)=EDGEC4) 

ENDIF 

DO  504  1=1,7 
NSC  I ) =0 

IF  (EDGE(I).LT.EDGECI+1) ) THEN 

DSIGrEUALI ( EDGE ( I ) , EDGE ( 1 + 1 ) , TSAI,  NSC  I ) ) 

IF  C DEBUG ) WRITE  Cl,*)  CH AR C I CHAR C ' A ' ) +1 - 1 ) , 

+ EDGE  C I ) , EDGE  C 1 + 1 ) , DSIG 

UAL=UAL+DSIG 
ENDIF 

504  CONTINUE 

UAL  = UAL*ALPHA**3*P32*AM0/ C PIE*EP*S3) 

IF  C DEBUG ) WRITE  Cl, 31)  NS, UAL 
31  FORMAT  C IX, 715, E10. 3) 

RETURN 

END 
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STITL  TSAI  --  Evaluate  the  integrand  of  Tsai's  formula 
SUBROUTINE  TSAI < CTHK, SIG ) 

C Subroutine  to  calculate  the  integrand  of  the  formula  for 

C the  brems.  cross  section  given  by  Tsai  in  SLAC-PUB-B49  (Jan 

C 1971)  (Formula  A. 24).  The  constants  in  the  commons  CONST 

C and  KIN  must  be  set.  Angles- in  radians#  energies  and 

C momenta  in  fm-1 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /CONST/  EM# EM2# AM# HBARC# ALPH# PIE# ALPHA, RTD# DTR 

COMMON  /KIN/  EEX#  TH, AMO#  S3#  P3#  CTH#  STH#  S32#  ES#  P32#  EP#  AMF, SDP3. SDP# 

+ U32#U3#U0#  U2#  CTHS#  CTHP#  STHP#  XMU 
COMMON  /CORREC/IWTILD#  WTIL1#  WTIL2#  ISOFT#  S0FT1 
C sin  of  photon  angle 

STHK  rDSQRT ( 1 . D0-CTHK**2) 

Energy  of  the  real  photon  (avoid  the  cancelation  of  mass  of  target 
squared ) 

W = ( EM2-SDP+AM0* ( ES-EP-XMU ) ) / ( U0-U3*CTHK ) 

C Momentum  part  of  s-p-k 

Q32=U32+W**2-2. D0*U3*W*CTHK 
Q3=DSQRT(Q32) 

C Length  of  q=s-p-k  (4-momentum  transfer  to  nucleus) 

Q2: (ES-EP-W)**2-Q32 
C 

A=W*(EP-P3*CTHP*CTHK) 

AA=W*(ES-S3*CTHS*CTHK ) 

B2= ( W*P3*STHP*STHK )**2 
X=DSQRT(A**2-B2) 

YrDSQRT (AA**2-B2) 

C The  constant  NU*(l/x-l/y)  of  Tsai 

XYNU=(AA+A)/(X*Y*(X+Y) ) 

C Calculate  the  integrand 

SIG=W*( F0RM2C  Q2)*C -A*EM2*(2. DO*ES* ( EP  + W ) + . 5D0*Q2 ) /X**3  - 
+ AA*EM2*(2. D0*EP*(ES-W)+. 5D0*Q2 )/Y#*3  - 2. DO  + 

+ 2.  D0*XYNU*(  EM2*  (SDP-W*>f<2  ) + SDP*  ( 2 . D0*ES*EP-SDP  + U*  ( ES-EP  ) ) ) + 

+ (2. D0*(ES*EP+ES*W+EP**2)+. 5D0*Q2-SDP-EM2)/X  - 
+ (2. D0*(ES*EP-EP*W+ES**2)+.5D0*Q2-SDP-EM2)/Y)  + 

+ F0RM1 (Q2)*( ( A/X**3+AA/Y**3 ) *EM2* ( 2. D0*EM2+02)  + 4 . DO  + 

+ 4 . DO*XYNU*SDP* ( SDP  — 2 . D0*EM2  ) + 

+ ( 1 . DO/X-1 . D0/Y)*(2. D0*SDP+2. D0*EM2-Q2) ) ) / ( 02*Q2* ( U0-U3*CTHK ) ) 

IF  (ISOFT. EQ.O)  RETURN 

SIG=SIG*S0FT1**( ALPHA* ( DLOG ( -Q2/EM2 ) -1 . DO) /PIE) 

RETURN 

END 


68 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 


•$T I TL  PWASPL  — Piece-wise  analytic  integration  of  Tsai 
SUBROUTINE  PWASPL(UAL) 

C We  assume  that  the  COMMON  KIN  has  been  setup  already 

C by  a call  to  INIKIN  and  that  the  common  PWAFF  has  been 

C initialized  with  a call  to  PWAFFI. 

‘ C Size  of  arrays  in  PWAFF 

PARAMETER  (MNPT=200) 

IMPLICIT  DOUBLE  PRECISION  (A-H/O-Z) 

DOUBLE  PRECISION  X ( MNPT ) , A ( MNPT ) , B ( MNPT ) , C ( MNPT ) , D ( MNPT ) 

DOUBLE  PRECISION  AP ( MNPT ), BP ( MNPT ), CP ( MNPT ), DP ( MNPT ) 

COMMON  /PWAFF / NPT , X, A , B , C , D , AP , BP , CP/ DP 

COMMON  /KIN/  EEX,TH,AM0,P1,P2, CTH, STH,P12,El,P22, E2, AMF/SDP3, SDP, 
+ U32#  U3»U0/U2»CTHS/CTHP  / STHP, XMU 
COMMON  /CONST/  EM,  EM2> AM/ HBARC,  ALPH, PIE/ ALPHA, RTD, DTR 
UAL=0. DO 

CALL  PWAINICXl/XMAX) 

CALL  ST ARTT ( 3 ) 

DO  30  1=1, NPT- 1 

IF  (Xl.LT.XCI) ) GO  TO  40 
30  CONTINUE 

I =NPT 

40  CALL  ST0PTC3) 

50  X2  = DMIN1 (X(I ) , XMAX) 

CALL  PWAEUACX1, X2, ACI ) , B( I ) , C< I >, D( I ) , 

+ AP< I ) , BP( IJ , CP( I)  , DP( I ) , XINT) 

UAL  =UAL  + X I NT 
1 = 1 + 1 
XI  =X2 

IF  (X2.LT.XMAX)  GO  TO  50 

UAL  = UAL*ALPHA*>K3*P22/(  PIE*E2*P1  ) 

RETURN 

END 
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$TITL  PWAINI*  PWAEUA  — Analytic  Tsai*  cubic  form  factor 
SUBROUTINE  PWAI N I < XM I N , XMAX ) 

PWAINI  should  be  called  once*  after  the  common  KIN  has 
been  initialized  in  order  to  calculate  constants  which 
do  not  change  when  only  THETA-k  Cor  q-MU)  changes. 

IMPLICIT  DOUBLE  PRECISION  CA-Z) 

COMMON  /PWAC0M/K1M2, K1M1, K10*K1P1*  K1P2*  K1P3*  K2M2*  K2M1, K20* 

+ K2P1, K2P2.K2P3, J1M2*J1M1*J10,J1P1*J1P2*  J2M2,  J2M1 * J20*  J2P1* 

+ J2P2*  L » LAMBDA 

COMMON  /CONST/  EM , EM2 * AM  * HBARC . ALPH * P I E * ALPHA  * RTD * DTR 
COMMON  /KIN/  EEX*TH, AMO*  PI*  P2*  CTH* STH, PI 2*  E 1 * P22 * E2*  AMF*SDP3, SDP* 
+ U32 * U3*  UO*  U2, CTHS » CTHP,  STHP*XMU 
LAMBDA=SDP-EM2 
KAPPA: E1-E2-XMU-LAMBDA/AM0 

CALL  PWACOEC LAMBDA*  KAPPA, P 1 , El  * P2*  E2* 1 . DO*  Al*  B1 * Cl*  ETA 1 * 

+ XMI N*  BET0/BET1* BET2*  TAU1*RH01*  BETOP*  BET  IP, BET2P ) 

CALL  PWACOE (LAMBDA, KAPPA*  P2* E2 * P 1 * E 1 * - 1 . DO * A2 * B2 * C2 * ETA2* 

+ XMAX*  DELO*  DELI  * DEL2* TAU2,RH02, DELOP*  DELIP,  DEL2P ) 

GAMMA: ( ( L AMBD A + EM2 ) * ( 2 . DO*E 1 *E2-LAMBDA ) +LAMBDA*K APPA**2+ 

+ ( LAMBDA  + EM2 ) *<  LAMBDA/ AMO  + XMU ) *KAPPA ) /LAMBDA**2 
ETA:1 . DO 

TAUP:-2. D0*(LAMBDA**2-2. D0*EM2* ( LAMBDA+EM2 ) ) /LAMBDA 
GAMMAP:2. D0*(LAMBDA**2-EM2*#2)/LAMBDA#*2 
ETAP :-2 . DO 

RHOP : C LAMBDA **2-2 . D0*EM2**2 )/LAMBDA**2 
RETURN 
C 

C PWAEUA  calculates  the  integral  of  Tsai  from  XI  to  XF  assuming 

C that  the  longitudinal  form  factor  (W2)  is  a cubic  defined 

C by  A*  B*  C*  and  D and  that  the  transverse  form  factor  (Wl) 

C is  a cubic  given  by  AP,  BP*  CP*  and  DP. 

ENTRY  PWAEUA (XI , XF , A , B * C * D , AP , BP * CP*  DP * UAL) 

CALL  PW AK I N ( X I * XF, LAMBDA*  KAPPA*  P1*E1*P2*E2*ETA1*A1*B1*C1* 

+ Dll, DIF, J10*  J1M1, J1P1*J1M2*  J1P2*K10*K1M1*K1P1*K1 M2 , K1P2, K1P3) 
CALL  PW AK I N ( X I * XF*  LAMBDA*  KAPPA*  P2*  E2*  PI*  El  * ETA2* A2*  B2*  C2* 

+ D2 1 * D2F i J20* J2M1, J2P1 , J2M2* J2P2*  K20 * K2M 1 * K2P 1 , K2M2* K2P2 , K2P3) 
BETA:2. D0*KAPPA*DSQRT( LAMBDA* (LAMBDA+2. D0*EM2) ) 

ZETA1I: (2. 0D0*LAMBDA*A1-B1 )*XI-2. D0*LAMBDA*B1 +C1 
ZETA1F: (2. 0D0*LAMBDA*A1-B1 )*XF-2. DO*LAMBDA*B 1 +C 1 
ZETA2I: (2. 0D0*LAMBDA*A2-B2)*XI-2. DO*L AMBDA*B2+C2 
ZETA2F: (2. 0D0*LAMBDA*A2-B2) *XF-2 . D0*LAMBDA*B2+C2 
C This  is  the  formula  as  originally  written  in  Maximon's  notes 

C L:-DLOG( ( ZET A2F+BETA*D2F ) / ( ZETA 1 F+BETA*D IF ) * 

C + (ZETA1I+BETA*D1I)/(ZETA2I+BETA*D2I) )/BETA 

C The  following  version  of  L is  rewritten  to  avoid  cancelation 

S 1 1 : 1 . DO 
EPS1 1:0. DO 

IF  (ZETA1I.LT. 0. DO ) THEN 

5 1 I : — 1 .DO 
EPS  1 1 :2 . DO 

END  I F 
S2I : 1 . DO 
EPS2 I :0 . DO 

IF  (ZETA2I .LT.O. DO)  THEN 

52  I :- 1 . DO 
EPS2I : 2 . DO 

END  I F 
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sir ; i . do 

EPS  1 F = 0 . DO 

IF  ( ZETA1F . LT . 0 . DO ) THEN 
S1F=-1 .DO 
EPS1F  = 2. DO 
ENDIF 
S2F=1 . DO 
EPS2F  = 0 . DO 

IF  CZETA2F.LT. 0. DO)  THEN 
S2F=-1 .DO 
EPS2F  = 2 . DO 
ENDIF 

DEN=2. D0*EM*P1*P2*KAPPA*STH 

L = C S IF  *DLOG ( (DABSCZETA1F )+BETA*DlF )/DEN)- 
+ SI I *DLOG ( ( DABS  C ZETA1 1 ) +BETA*D1 1 ) /DEN ) - 
+ S2F  *DLOG ( C DABS CZETA2F)+BETA*D2F)/DEN)+ 

+ S2I*DL0GC ( DABS  C ZET A2 1 > + BETA*D2I ) /DEN ) + 

+ EPS IF  #DLOG ( DABS  C XF-2 . DO*LAMBDA) )- 
+ EPS1 I*DLOGC DABSCXI-2. DO*LAMBDA) )- 
+ EPS2F*DL0G( DABSCXF-2 . DO^LAMBDA) )+ 

+ EPS2I*DL0G(DABS(XI-2.D0*LAMBDA)) ) /BET  A 
CALL  PWAINTCXI, XF, BETO, BET1 , BET2,  DELO/ DELI / DEL 2/ ETA/ 

+ TAU1/ T AU2 / RH01 / RH02/  GAMMA/ A/ B/ C/  D/ I ) 

CALL  PWAINTCXI, XF, BETOP, BET IP, BET2P, DELOP,  DELIP, DEL2P, ETAP, 
+ TAUP, TAUP/RHOP, RHOP, GAMMAP, AP, BP, CP, DP, IP) 

OAL=I+IP 

RETURN 

END 
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$TITL  PWACOE  — Initialize  constants  constant  with  q-MU 

SUBROUTINE  PWACOE ( LAMBDA , KAPPA , PA , EA , PB , EB , SON , A , B, C » ETA, 

+ X,  BETO , BET1 * BET  2, T AU , RH0,BET0P,BET1P,BET2P) 

C Note  that  the  argument  SGN  will  be  either  +1  or  -1.  This 

C permits  one  to  use  the  same  formulae  to  calculate  two  sets 

C of  variables  necessary  for  the  PWA  calculation. 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

COMMON  /CONST / EM, EM2, AM/ HBARC, ALPH, PIE, ALPHA, RTD, DTR 
COMMON  /KIN/  EEX, TH, AMO, PI, P2, CTH, STH, PI 2, El, P22, E2, AMF , SDP3, SDP, 
+ U32, U3, U0,U2, CTHS, CTHP, STHP,XMU 
Ar (PB+SGN*(EA*PB-EB*PA*CTH)/AM0)**2+(EM*PA*STH/AM0)**2  ' 

B=2. D0*PA*(PA*PB-(EA*EB-EM2)*CTH)*(PB+SGN*(EA*PB-EB*PA*CTH)/ 

+ AMO ) +SGN*2 . D0*EM2*PA**2*(EA-EB)*STH**2/AM0+2. DO*XMU* 

+ ( -SGN*PB*  (EA*P8-EB*PA*CTH)  -LAMBDA*  (LAMBDA+2.D0*EM2)  / AMO ) 

Cs4 . DO* ( PA*LAMBDA-SGN*XMU* ( PA*EB-EA*PB*CTH) )**2  + 

+ (2. D0*EM*XMU*PB*STH)**2 
ETA= ( (2. D0*EM*PA*PB*KAPPA*STH)/A)**2 
X=2. DO* (LAMBDA* ( EA-EB-SGN*XMU ) * ( E A-EB-U3 ) )/ 

+ (1 . DO+SGN*(EA-EB-U3)/AMO) 

PAB=SGN*PA*(PA-PB*CTH) 

EAMU=EA*(EA-SGN*XMU) 

EAAMOrl . D0+SGN*2 . DO*E A/AMO 

BETO =8. D0*EM2*EAMU*LAMBDA*( PAB- ( EA+EB )*XMU) 

BET  1 z-EM2* ( 4 . DO*EAMU* C PAB- ( EA+EB ) * ( KAPPA  + XMU ) ) + 

+ 2. DO*L AMBDA* CPAB-C EA+EB )*XMU)*EA AMO) 

BET2=EM2*CPAB-( EA+EB)* (KAPPA+XMU) ) *EAAMO 
TAU= (LAMBDA**2-4 . D0*EM2*EA*EB-SGN*2 . DO*XMU* 

+ (EM2*(EA-EB-SGN*XMU)+EB*LAMBDA) ) /LAMBDA 
RHO=  ( (LAMBDA+EM2 )*(2. DO*EA*EB-L AMBDA ) +LAMBDA*KAPPA**2+ 

+ (LAMBDA+EM2)*(LAMBDA/AM0+XMU)*KAPPA)/LAMBDA**2+ 

+ EM2/AM0**2+. 5D0+SGN*E A/AMO 
BET0P=-8. D0*EM2**2*LAMBDA* (PAB- (EA+EB )*XMU) 

BET  IP: 4 . D0*EM2* ( EM2* ( PAB- ( EA  + EB )*( KAPPA  + XMU ) ) + 

+ LAMBDA*(PAB-(EA+EB)*XMU) ) 

BET2P: -2 . D0*EM2*( PAB- ( EA+EB )*( KAPPA+XMU) ) 

RETURN 

END 
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$TITL  PWAKIN  — Evaluate  constants  which  change  with  q-MU 

SUBROUTINE  PW AK I N ( X I , XF , LAMBDA , KAPPA , PA , EA , PB , EB , ETA, A, B, C, 

+ DI , DF, JO, JM1, JP1, JM2, JP2,K0, KM1,  KP1, KM2,  KP2,  KP3 ) 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

COMMON  /CONST/  EM, EM2 , AM , HB ARC , ALPH , P IE , ALPHA , RTD , DTR 
COMMON  /KIN/  EEX, TH, AMO, PI, P2, CTH, STH, PI 2, El , P22, E2, AMF , SDP3 , SDP , 
+ U32»U3,U0»  U2 , CTHS , CTHP, STHP » XMU 
CHI=B/A 

D I : DSQRT ( A* ( ( XI -CH I ) **2+ETA ) ) 

DF  r DSQRT ( A* ( ( XF-CH I )**2+ETA) ) 

SQTAr DSQRT ( A ) 

SQTCzDSQRT(C) 

C The  following  is  the  formula  as  originally  written 

C JO=DLOG( ( A*XF-B+DF*SQTA)/(A*XI-B+DI#SQTA) )/SQTA 

C This  is  the  version  as  corrected  to  reduce  "cancelation" 

SI  = 1..  DO 

IF  (A*XI-B.LT.O. DO)  SI=-1.D0 
SF;1 .DO 

IF  ( A*XF-B.LT.O. DO)  SF=-1.D0 
DEN=2. DO*EM*PA*PB*KAPPA*STH 
J0=  (SF*DLOG(  ( DABS  ( A*XF-B)+DF *SQTA  ) /DEN  ) - 
+ SI*DLOG( (DABS( A*XI-B)+DI*SQTA)/BEN) )/SQTA 
C The  following  is  the  formula  as  originally  written 

C JM1 =-DLOGCXI/XF*(-B*XF+C+SQTC*DF)/ 

C + (-8*XI+C+SQTC*DI) )/SQTC 

C This  version  is  corrected  to  avoid  cancelation 

SI=1.D0 

IF  (B*XI-C.LT.O. DO)  SI=-1.D0 
SF= 1 . DO 

IF  (B*XF-C.LT.O.DO)  SF=-1.D0 

JM1 z ( SF *DLOG ( CDABSC  B*XF-C)+SQTC*DF)/(XF*DEN)  )- 
+ SI *DLOG ( ( DABS( B*XI-C ) +SQTC*DI )/<XI*DEN) ) )/SQTC 
JP1= ( DF-DI ) / A + CH I * J 0 
JM2= (-DF/XF+DI/XI+B*JM1 )/C 

JP2 - (XF*DF-XI*DI+3. DO*B* JP1-C* JO ) / ( 2 . DO* A ) 

KOr  C CA*XF-B)/DF-(A*XI-B)/DI )/(ETA*A**2) 

KPl=-(  (-B>i<XF  + C)/DF-(-B*XI+C)/DI  ) / ( ETA*A**2  ) 

KP2= ( JO+2. D0*B*KP1-C*K0)/A 

KP3:  (XF*  + 2/DF-XI**2/DI+3.  D0*B*KP2-2.  D0*C>*<KP1  ) /A 
KM1= ( 1 . DO/DF-1 . DO/D 1+ JM1 +B*KO)/C 

KM2= (-1 . D0/(XF*DF)  + 1 . DO/(XI*DI )+3. D0*B*KMl-2.  DO*A*KO)/C 

RETURN 

END 
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STITL  PWA1NT  — Determine  integral  of  long,  or  trans  term. 

SUBROUTINE  PWAINT ( XI , XF, BETO, BET1 , BET2, DELO, DELI , DEL2, ETA, 

+ T AU1 »TAU2,RH01,RH02, GAMMA, A, B,  C,  D, I ) 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

COMMON  /PWAC0M/K1M2,K1M1,  K10,K1P1,  K1P2,  K1P3, K2M2 , K2M1 , K20 , 

+ K2P1 , K2P2, K2P3, J1M2, J1M1 , J 1 0, J1P1#J1P2, J2M2 » J2M 1 * J20, J2P 1 , . 

+ J2P2, L, LAMBDA 

COMMON  /KIN/  EEX, TH, AMO, PI , P2, CTH, STH, P12, El, P22, E2, AMF , SDP3, SDP, 
+ U32,U3,U0,U2, CTHS , CTHP,  STHP,XMU 
CALL  PWAE INC  BETO, BET1 , BET2, K1M2, K1M1 ,K10,K1P1,K1P2,K1P3, 

+ J1 M2, J1M1,J10,J1P1,J1P2, A, B,C,D, 111,131,141) 

CALL  PWAE INC  DEL 0, DELI , DEL2, K2M2, K2M1 , K20, K2P1 , K2P2, K2P3, 

+ J2M2, J2M1 , J20, J2P1 , J2P2, A, B,C, D, 112, 132, 142) 

I2=A*< 1 . DO/XI-1 . DO/XF)+B*DLOG(XF/XI )+C*(XF-XI )+ 

+ . 5D0*D*(XF**2-XI**2) 

15: ( A+B*2. D0*LAMBDA+C*4. D0*LAMBDA**2+D*8 . D0*LAMBDA**3 ) *L- 
+ ( B + 2 . D0*LAMBDA*C+4.D0#LAMBDA*#2#D)*(J10-J20)- 
+ CC+2. DO*LAMBDA*D)*( J1P1-J2P1 )- 
+ D*( J1P2-J2P2) 

I:-I11-I12-ETA#I2/U3+TAU1*I31-TAU2*I32- 
+ RH01*I41+RH02*I42-GAMMA*I5 
RETURN 
END 
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1 STITL  PWAEIN  --  Evaluate  some  terms  in  the  analytic  integral 

2 SUBROUTINE  PWAEINC BETO, BET1 , BET2, KM2, KM1 / KO, KP1 , KP2, KP3» 

3 + JM2, JM1. JO, JP1, JP2, A, B, C, D, 11/ 13, 14 ) 

4 IMPLICIT  DOUBLE  PRECISION  (A-Z) 

5 Il=BET0*A*KM2+( BETO*B+BET 1 *A )*KM1+ 

6 + (BET0*C+BET1*B+BET2*A)*K0+(BET0*D+BET1*C+BET2#B)*KP1+ 

7 + <BET1*D+BET2*C)*KP2+BET2*D*KP3 

8 I3=A*JM2+B*JM1+C*J0+D*JP1 

9 I4=A*JM1+B*J0+C*JP1+D*JP2 

10  RETURN 

11  END 
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ST1TL  SPCE  — Series  expansion  of  Spence  function 
FUNCTION  SPCE ( X ) 


Series  expansion  of  the  Spence  function  for  range  0 
Note  that  the  number  of  terms  used  in  the  series  is 
to  auo  i d floating  underflow. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z> 

DOUBLE  PRECISION  SPCE 
SPCEsO. ODO 

IF  (X.LE.O.DO)  RETURN 
TOP  = 1 . ODO 

N=1  . DO-25. DO/DLOGCX) 

DO  10  1=1, N 

TOP=TOP*X 

SPCE=SPCE+T0P/I**2 
10  CONTINUE 
RETURN 
END 


to  0.5 
chosen 
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1 STITL  PWAEIN  — Evaluate  some  terms  in  the  analytic  integral 

2 SUBROUTINE  PWAEINC BETO, BET1 , BET2, KM2, KM1 , KO, KP1 , KP2, KP3, 

3 + JM2»  JM1*J0>JP1*JP2/A»B»C,D» 1 1 > I3»  14) 

4 IMPLICIT  DOUBLE  PRECISION  CA-Z) 

5 Il=BETO*A*KM2+(BETO*B+BETl*A)*KMl+ 

6 + <BET0*C+BET1*B+BET2*A)*K0+(BET0*D+BET1*C+BET2*B)*KP1+ 

7 + <BET1*D+BET2*C)*KP2+BET2*D*KP3 

8 I3rA*JM2+B*JMl+C*J0+D*JPl 

9 I4=A*JM1+B*J0+C*JP14D*JP2 

10  RETURN 

11  END 
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SPCE  --  Series  expansion  of  Spence  function 
FUNCTION  SPCE  C X ) 

Series  expansion  of  the  Spence  function  for  range  0 
Note  that  the  number  of  terms  used  in  the  series  is 
to  avoid  floating  underflow. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  SPCE 
SPCE=0. ODO 

IF  (X.LE.O.DO)  RETURN 
TOP  r 1 . ODO 

N=1 . DO-25. DO/DLOGCX) 

DO  10  1=1, N 

TOP=TOP*X 

SPCE=SPCE+T0P/I**2 

CONTINUE 

RETURN 

END 


to  0.5 
chosen 


1 $TITL  SPENCE  — The  Spence  function 

2 FUNCTION  SPENCE  C X ) 

3 C Program  for  the  computation  of  the  Spence  function 

4 C values  of  the  argument  x 

5 IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

6 DOUBLE  PRECISION  SPENCE 

7 DATA  PI23, PI26/3. 289868134D0, 1 . 644934067D0/ 


a 

IF 

(X.GT.2.0D0) 

GO 

TO 

1 

9 

IF 

(X. GT. 1 .ODO) 

GO 

TO 

2 

10 

IF 

(X.GT.0.5D0) 

GO 

TO 

3 

n 

IF 

<X. LT. -1 . DO) 

GO 

TO 

4 

12 

IF 

(X.LT.O.ODO) 

GO 

TO 

5 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 


SPENCE=SPCE<X) 

RETURN 

Y = 1 . ODO/X 

SPENCE  = PI 23-SPCE < Y ) - . 5D0*DL0G ( X ) **2 
RETURN 

Y=(X-1 . ODO)/X 
Yl:  (X-l  ..  ODO  ) * Y 

SPENCE=PI26+SPCE(Y)-0. 5D0*DL0G  <X)*DL0G(Y1) 

RETURN 

Y=1 .ODO-X 

SPENCE  = PI 26-SPCE ( Y ) -DLOG  < X ) *DLOG( Y ) 

RETURN 
Y=1 . ODO-X 
Yl=l .ODO/Y 

SPENCE =-PI 26+SPCE< Yl )-0. 5D0*DL0G < Y ) *DLOG < X*X/Y ) 
RETURN 

Y zX/ (X— 1 . ODO) 

Yl=-X+1 . ODO 

SPENCE  = -SPCE  C Y ) -0 . 5 DO* DLOG ( Yl ) 

RETURN 

END 


for  all  real 
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*TITL  EUALI  — Simpson's  rule  integration 
FUNCTION  EUALI ( A / B / FUNC  / NS) 

C Function  to  evaluate  the  integral  of  the  function  FUNC 

C betuieen  6 and  B.  The  number  of  steps  starts  at  NSTEP 

C as  defined  in  EUALI  and  then  goes  up  by  a factor  of  two 

C until  the  integral  stops  changing  by  more  than  a factor 

C of  EPS.  In  this  version/  which  uses  Simpson's  rule/  the 

C new  values  needed  when  NSTEP  increases  are  not  all  re- 

C calculated.  Ualues  computed  during  the  calculation  of 

C the  old  (inaccurate)  value  of  the  integral  are  re-used. 

IMPLICIT  DOUBLE  PRECISION  (A-H/O-Z) 

EXTERNAL  FUNC 
DOUBLE  PRECISION  EUALI 
COMMON  /EU I COM/  NSTEP/ EPS 
CALL  FUNC( A/ EUALA) 

CALL  FUNC ( B/ EUALB ) 

NSrNSTEP-1 
STEPO=  ( B- A ) /"NSTEP 
EUAL2=0. DO 
DO  900  1=1 /NS 

CALL  FUNC (A+I*STEPO/ U) 

EUAL2  = EUAL2  + U 

900  CONTINUE 
NS=NSTEP 

C 

C Start  of  loop 

4 STEP: . 5D0*STEP0 

EUAL4=0. DO 
DO  901  1=1/ NS 

CALL  FUNC (A+STEP+C 1-1 ) *STEPO/ U ) 

EUAL4=EUAL4+U 

901  CONTINUE 

EU=  (EUALA  + 2.  D0*EUAL2  + 4 . D0*EUAL4+EUALB  ) *STEP0/'6  . DO 
IF  (NS. NE. NSTEP)  THEN 

IF  (EPS. LE. 0. DO)  THEN 

IF  (DABS(EU-EUALIO) . LE . DABS ( EU#EPS ) ) GO  TO  6 
CD  WRITE  (1,*)  A/ B/ DABS (EU-EUALIO) / DABS( EU*EPS ) 

ELSE 

IF  ( D ABS ( EU-EUAL 1 0 ) . LE . EPS ) GO  TO  6 
CD  WRITE  (1,*)  A/ B/ DABS ( EU-EUAL 1 0 ) / EPS 

END  IF 
END  I F 
. EUALIOrEU 
EUAL2  = EUAL4  + E VAL2 
NS  = 2*NS 
STEPO=  STEP 
GO  TO  4 

6 EUAL I =EU 

RETURN 
END 
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JTITL  EUALI  --  Gaussian  integration 
FUNCTION  EUALI(A,B,FUNC,NS) 

C Function  to  evaluate  the  integral  of  the  function  FUNC 

C between  A and  B.  The  number  of  intervals  starts  at  NSTEP 

C as  defined  in  EUICOM  and  then  goes  up  by  a factor  of  two 

C until  the  integral  stops  changing  by  more  than  a factor 

C of  EPS.  In  this  version,  NGAUSS-point  Gaussian  integration 

C is  used  (for  particularly  tough  functions). 

PARAMETER  (NGAUSSr 16) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

EXTERNAL  FUNC 

DOUBLE  PRECISION  EUALI , 2 ( NGAUSS ) > W ( NGAUSS ) 

COMMON  /EUICOM/  NSTEP, EPS 
NS=NSTEP 
C 

C Start  of  loop 

4 STEP  = ( B— A ) /NS 

EUrO.DO 
DO  902  1=1, NS 

AX=A+( 1-1 )*STEP 
BX= AX+STEP 

CALL  GSET  ( AX, BX, NGAUSS , 2, W ) 

DO  903  J=l, NGAUSS 

CALL  FUNC (2(J),U) 

EU=EU+U*W( J) 

903  CONTINUE 

902  CONTINUE 

IF  (NS. NE. NSTEP)  THEN 

IF  (EPS. LE. 0. DO)  THEN 

IF  ( DABS ( EU-EUALIO ) . LE . DABS ( EU*EPS ) ) GO  TO  6 
CD  WRITE  (1,*)  A, B, DABS(EU-EUALIO) , DABS(EU*EPS) 

ELSE 

IF  (DABS(EU-EUALIO) . LE . EPS)  GO  TO  6 
CD  WRITE  (1,*)  A, B, DABS ( EU-EUAL I 0 ) , EPS 

ENDIF 
END  IF 
EUALIOrEU 
NS=  2*NS 
GO  TO  4 

6 EUALI=EU 

RETURN 
END 
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1 

ST1TL 

GSET 

— N-point  Gauss  zeros  and  weights 

z 

SUBROUTINE  GSET  ( AX , BX , NX , Z , W ) 

3 

C 

N-point  Gauss  zeros  and  weights  for  the 

4 

c 

stored  in  arrays  Z and  W respectively. 

5 

c 

6 

c 

This 

version  converted  from  the  version 

7 

c 

for  the  CDC-7600  to  a double  precision  v 

8 

c 

mach ine . 

9 

c 

10 

IMPLICIT  DOUBLE  PRECISION  CA-H,0-Z> 

11 

DOUBLE  PRECISION  Z(21»W(21>A(2731/X(2731 

12 

INTEGER  KTABC96) 

13 

c 

14 

c 

Table 

of  initial  subscripts  for  N=2(l)16 

15 

DATA 

KTAEK21/1/ 

16 

DATA 

KTABOl/2/ 

17 

DATA 

KTAB  C 4 ) /4/ 

18 

DATA 

KTAB  ( 5 ) /6/ 

19 

DATA 

KT AB  C 6 ) /9/ 

20 

DATA 

KTAB(71/12/ 

21 

DATA 

KTAB ( 8 1 / 16/ 

22 

DATA 

KTAB  ( 9 ) /20/ 

23 

DATA 

KTAB( 101/25/ 

24 

DATA 

KTAB( 11 1/30/ 

25 

DATA 

KTAB( 121/36/ 

26 

DATA 

KTAB( 131/42/ 

27 

DATA 

KTAB(141/49/ 

28 

DATA 

KTAB ( 151/56/ 

29 

DATA 

KTAB  C 161/64/ 

30 

DATA 

KTAB ( 20 ) /72/ 

31 

DATA 

KTAB ( 24 ) /82/ 

32 

DATA 

KTAB ( 28 1 /B2/ 

33 

DATA 

KTABC321/94/ 

34 

DATA 

KTAB ( 36 ) /94/ 

35 

DATA 

KTAB (401/110/ 

36 

DATA 

KTAB(441/110/ 

37 

DATA 

KTAB(481/130/ 

38 

DATA 

KT  AB ( 52 1/130/ 

39 

DATA 

KTAB(56)/130/ 

40 

DATA 

KTAB(601/130/ 

41 

DATA 

KTAB (64 1/154/ 

42 

DATA 

KTAB(681/154/ 

43 

DATA 

KTABC721/154/ 

44 

DATA 

KTAB(761/154/ 

45 

DATA 

KTAB(801/ie6/ 

46 

DATA 

KTAB(841/186/ 

47 

DATA 

KTABC881/186/ 

48 

DATA 

KTABC921/186/ 

49 

DATA 

KT  AB ( 96 1/226/ 

50 

c 

51 

c 

Table 

■ of  abscissae  (XI  and  weights  (A!  fi 

52 

c 

53 

c 

N = 2 

54 

DATA 

X( 1 1 /O. 577350269189626D0  /,  A(ll/1 

55 

c 

N = 3 

56 

DATA 

X(21/0. 774596669241483D0  /,  A(21/0 

57 

DATA 

X(31/0. OOOOOOOOOOOOOOODO  /,  A(3)/0 

A( 31/0. 888888888888B89D0 
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58 

C 

N = 4 . 

59 

DATA 

X 

<4)/0. 861 13631 1594053D0 

A ( 4 ) SO . 

60 

DATA 

X 

(5)/0. 3399810435B4856D0 

A ( 5 ) SO . 

61 

C 

N = 5 

62 

DATA 

X 

( 6 ) /"0 . 9061 798 45938 664 DO 

/ , 

A ( 6 ) so . 

63 

DATA 

X 

<7)/0. 5384693101 05683D0 

/ , 

A ( 7 ) so . 

64 

DATA 

X 

( 8 ) / O . 000 000 000 000 000 DO 

/ , 

A ( 8 ) SO . 

65 

C 

N = 6 

66 

DATA 

X 

( 9) so. 9324695 1 4203 152D0 

/ , 

A ( 9 ) SO . 

67 

DATA 

X 

( 10)/0. 661209386466265D0 

/ , 

A ( 1 0 ) SO 

68 

DATA 

X 

(11  )/0. 23861 9186083 197D0 

/ , 

A(  11  )sO 

69 

C 

N = 7 

70 

DATA 

X 

( 12)/0. 949107912342759D0 

/ , 

A ( 1 2 ) SO 

71 

DATA 

X 

< 13 )S0. 741531 185599394 DO 

/ , 

A ( 1 3 ) /O 

72 

DATA 

X 

( 14 )/0. 405845 151 377397D0 

/ , 

A ( 141/0 

73 

DATA 

X 

( 15 ) SO . 000000000 000000 DO 

/ , 

A( 15)/0 

74 

C 

N = 8 

75 

DATA 

X 

( 1 6 ) /O . 96 028 9 85 6 49 7 53 6 DO 

/ , 

A( 16)/0 

76 

DATA 

X 

( 17)/0. 7966664774 13627D0 

/, 

A ( 17 ) SO 

77 

DATA 

X( 18) SO. 52553240991 6329D0 

/ , 

A ( 18) SO 

78 

DATA 

XC19)/0. 183434642495650D0 

/ , 

A ( 1 9 ) SO 

79 

C 

N = 9 

80 

DATA 

X C 20 ) SO . 9681 60239507 626D0 

/ , 

A ( 20 ) SO 

81 

DATA 

X 

(21 1/0. 836031 107326636D0 

/ , 

A ( 2 1 ) SO 

82 

DATA 

X 

( 22 ) /O . 613371 4 3 270 0590 DO 

s , 

A ( 22 ) SO 

83 

DATA 

X 

( 23 ) /O . 324 25 34 2 3403 809 DO 

/ , 

A ( 23) SO 

84 

DATA 

X 

( 24 ) /Q . OOOOOOOOOOOOOOODO 

st 

A ( 24 ) SO 

85 

C 

N=10 

86 

DATA 

X 

(25)/0. 97390652851 7 172D0 

s , 

A(25)/0 

87 

DATA 

X 

(26)/0. 8650633566S8985D0 

St 

A ( 26 ) SO 

88 

DATA 

X 

(27)/0. 679409568299024 DO 

St 

A ( 27 ) SO 

89 

DATA 

X 

(28)/0. 433395394 129247 DO 

s , 

A ( 28 ) SO 

90 

DATA 

X 

(29)/0. 148874338981631D0 

Sr 

A(29)/0 

91 

c 

N = 1 1 

92 

DATA 

X 

( 30 )/0. 978228658 146057 DO 

S , 

A ( 30 ) SO 

93 

DATA 

X 

(31 )/0. 88706259976S095D0 

St 

A ( 31 )S0 

94 

DATA 

X 

(32)/0. 730152005574049D0 

St 

A ( 32 ) SO 

95 

DATA 

X 

(33)/0. 51 9096 12920681 2D0 

s , 

A( 33) SO 

96 

DATA 

X 

(34)/0. 269543155952345D0 

St 

A(34)/0 

97 

DATA 

X 

( 35 ) /O . OOOOOOOOOOOOOOODO 

s , 

A ( 35 ) SO 

98 

c 

N = 12 

99 

DATA 

X 

(36)/0. 981 5606342467 1 9D0 

s , 

A ( 36 ) SO 

100 

DATA 

X 

(37)/0. 90411 7 25 6 37 04 75 DO 

St 

A ( 37 ) SO 

101 

DATA 

X 

(38)/0. 7699026741 94 305D0 

St 

A ( 38 ) SO 

102 

DATA 

X 

(39)/0. 58731 79542866 1 7 DO 

s , 

A ( 39 ) SO 

103 

DATA 

X 

(40)/0. 36783 14989981 80 DO 

s , 

A ( 40 ) SO 

104 

DATA 

X 

(41 )/0. 12523340851 1469D0 

s , 

A ( 4 1 ) SO 

105 

c 

N = 13 

106 

DATA 

X 

( 42 ) /O . 984 18305471 8588D0 

s , 

A(42)/0 

107 

DATA 

X 

(43)/0. 91 7 598 399 22 2 97 8 DO 

St 

A ( 43 ) SO 

108 

DATA 

X 

(44 )/0. 80 1 57809073331 ODO 

s , 

A ( 44 ) SO 

109 

DATA 

X 

(45)/0. 642349339440340D0 

s , 

A ( 45 ) SO 

110 

DATA 

X 

(46)/0. 44 849275 1036447 DO 

s , 

A ( 46 ) SO 

111 

DATA 

X 

( 47 ) /O . 2304583 15955 135D0 

s , 

A ( 47 ) SO 

112 

DATA 

X ( 48 ) SO . 00 000 000 0000 000 DO 

s , 

A ( 48 ) SO 

113 

c 

N = 14 

114 

DATA 

X ( 49 )/0. 9862838086968 12D0 

s , 

A ( 49 ) SO 

, 101228536290376D0 
222381 034453374 DO 
313706645877887D0 
3626 8 3783 37 83 62 DO 

08 1 27438836 1574D0 


06667 1344308688D0 


2692667 19309996D0 


272925086777901  DO 


1 069393259953 18D0 
160078328543346D0 
203167426723066D0 
233492536538355D0 


✓ 

/ 

✓ 

/ 

✓ 

✓ 

/ 

/ 

/ 

/ 

✓ 

/ 

/ 

/ 

S 

s 

s 

s 

s 

s 

s 


040484004765316DO  / 
092121499837728D0  / 
138873510219787D0  / 
1781459B0761946D0  ✓ 
20781 6047536889D0  ✓ 
2262831 80262897 DO  ✓ 
23255 1553230874 DO  ✓ 

035119460331752D0  ✓ 
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115 

DATA 

X ( 50 ) /O . 928 4 34 08 366 3 57 4 DO 

/ , 

A(50)/0. 

116 

. 

DATA 

X ( 5 1 )/0. 827201 31 5069765D0 

/ , 

A( 51 ) /O. 

117 

DATA 

XC52)/0. 68729290481 1685D0 

/ , 

A(52)/0. 

118 

DATA 

X ( 53 )/0. 5 15248636358 154 DO 

/ , 

A(53)/0. 

119 

DATA 

X(54)/0. 3191 1 2 36 8927 890 DO 

/ , 

A(54)/0. 

120 

DATA 

X ( 55 ) /O . 108 05 4 94 8707 34 4 DO 

/ , 

A(55)/0. 

121 

C 

N = 1 5 

122 

DATA 

XC56)/0. 98799251 8020485D0 

/ , 

A(56)/0. 

123 

DATA 

X ( 57 ) /O . 937 27 3392 400706 DO 

/ , 

A(57)/0. 

124 

DATA 

X ( 58 )/0. 8482065834 10427D0 

/, 

A(58)/0. 

125 

DATA 

X(59)/0. 724417731360170D0 

/, 

A ( 59 ) /O . 

126 

DATA 

X(60)/0. 570972 17260B539D0 

// 

A( 60) /O. 

127 

DATA 

X ( 6 1 )/0. 394151347077563D0 

/ , 

AC61 )/0. 

128 

DATA 

X ( 62 ) /O . 201 194093997 435 DO 

/ , 

A(62)/0. 

129 

DATA 

X ( 63 ) /O . OOOOOOOOOOOOOOODO 

/ , 

A(63)/0. 

130 

C 

N=16 

131 

DATA 

X ( 64 ) /O . 98940093499 1650D0 

/ , 

A(64)/0. 

132 

DATA 

X ( 65  ) /O  . 94  4 57 5 02 307 323 3 DO 

/, 

A(65)/0. 

133 

DATA 

X(66)/0. 86563 1202387832D0 

// 

A ( 66 ) /O . 

134 

DATA 

X(67)/0. 755404408355003D0 

/, 

A ( 67 ) /O . 

135 

DATA 

X ( 68 )/0. 6 1787624440264 4D0 

/, 

A ( 68 ) /O . 

136 

DATA 

X ( 69 ) /O . 45801 677 7 657 22 7 DO 

/ , 

A ( 69 ) /O . 

137 

DATA 

X(70)/0. 281 60355077 9259D0 

/ , 

A<70)/0. 

138 

DATA 

X ( 7 1 >/0. 09501 2509837637D0 

/ , 

A ( 7 1 )/0. 

139 

c 

N:20 

140 

DATA 

X < 72 )/0. 993 128599 185094 DO 

/ , 

A ( 72 ) /O . 

141 

DATA 

X(73)/0. 963971 9272779 13D0 

/ , 

A ( 73 ) /O . 

142 

DATA 

X ( 74 ) /O . 912234428251 325D0 

/ , 

A(74)/0. 

143 

DATA 

X(75)/0. 839 11 6971 8222 18D0 

/ , 

A(75)/0. 

144 

DATA 

X ( 76 )/0. 74 633 1906460 150D0 

/, 

A ( 76 ) /O . 

145 

DATA 

XC77)/0. 6360536807265 15D0 

/, 

A ( 77 ) /O . 

146 

DATA 

X(79)/0. 51 0867001 950827D0 

/, 

A(78)/0. 

147 

DATA 

X(79)/0. 37 37060887 154 19D0 

// 

A(79)/0. 

148 

DATA 

X(80)/0. 227785851 14 1645D0 

/, 

AC  80/0. 

149 

DATA 

XC81 )/0. 076526521 133497D0 

// 

AC81 )/0. 

150 

c 

N r 24 

151 

DATA 

X ( 82 ) /O. 995 187 2 1999702 IDO 

/, 

A ( 82 ) /O . 

152 

DATA 

XC83)/0. 97472855597 1309D0 

// 

AC  83) /O. 

153 

DATA 

X(84)/0. 93827 4552002732D0 

/, 

A C 84 ) /O . 

154 

DATA 

X(85)/0. 8864 1552700440  IDO 

/, 

AC  85)/0. 

155 

DATA 

X(86)/0. 82000 1985973902D0 

/, 

A C 86 ) /O . 

156 

DATA 

X ( 87 )/0. 740124 191 578554 DO 

/, 

A C 87 ) /O . 

157 

DATA 

X(88)/0. 64809365 1936975D0 

// 

AC  88) /O. 

158 

DATA 

X(89)/0. 545421 47 1388839D0 

/, 

A C 89 ) /O . 

159 

DATA 

XC90)/0. 4 337 93507 62604 5D0 

/, 

A C 90 ) /O . 

160 

DATA 

X(91 )/0. 315042679696163D0 

/, 

AC91 )/0. 

161 

DATA 

X ( 92 ) /O . 19111 886747  36 1 6D0 

/, 

A C 92 ) /O . 

162 

DATA 

X(93)/0. 064056892862605D0 

/, 

AC93)/0. 

163 

c 

N = 32 

164 

DATA 

X(94)/0. 99726386 184 948  IDO 

/, 

AC94)/0. 

165 

DATA 

X(95)/0. 98561 151 1545268D0 

/, 

A C 95 ) /O . 

166 

DATA 

X(96)/0. 964762255587506D0 

/, 

A C 96 ) /O . 

167 

DATA 

X(97)/0. 934 906075937 739D0 

/, 

A C 97 ) /O . 

168 

DATA 

X(98)/0. 896321 155766052D0 

/, 

AC  98)/C. 

169 

DATA 

X ( 99 ) /O . 84936761 3 73 2 56 9 DO 

/ , 

AC99)/0. 

170 

DATA 

XC 1 00 )/0. 7 944 8379596794 2D0/, 

AC  100) /O 

171 

DATA 

X( 101 )/0. 7321821 18740289D0// 

AC  101 )/0 

000158007 159 7 60 DO  / 
121518570687903D0  / 
1572031 67 1581 94 DO  / 
1 85 5 38 397 4 77 9 38 DO  / 
205198463721296D0  / 
2152638534631 58 DO  / 

0307532419961 17 DO  / 
070366047488 108D0  / 
107 159220467 172D0  / 
1 39570677926 154D0  / 
1 662692058 16994 DO  / 
186161 0000 155 62 DO  / 
1984314853271 11D0  / 
20257824 1925561 DO  / 

027 152459411754D0  / 
062253523938648D0  / 
095158511682493D0  / 
124628971 255534 DO  / 
1 495959888 16577D0  / 
1691 565 1 9395003D0  / 
1 826034 15044924D0  / 
1 894506 10455069D0  / 

0176140071 39 1 52 DO  / 
040601 4 29 8 00 3 86 DO  / 
0626720483341 09 DO  / 
083276741 576704 DO  / 
1019301 19817240D0  / 
1 18194531961518D0  / 
131 6886384491 76D0  / 
142096109318382D0  / 
149172986472603D0  / 
1 52753387 130725D0  / 

01234 1229799987D0  / 
028531 3B8628933D0  / 
0442774388174 19D0  / 
059298584915436D0  / 
07334648141 1080D0  / 
086190161531953D0  / 
0976186521041 13D0  / 
1074442701 15965D0  / 
1 15505668053725D0  / 
121670472927803D0  X 
1 25837 4 56346828D0  / 
1 27930 195346752D0  / 

00701 8610009470D0  / 
016274394730905D0  / 
025 3 92 065 3 09 2 62 DO  / 
034273862913021D0  / 
0428 35 8980222 26 DO  / 
050998059262376DO  / 

. 058684  09  347 8 535 DO^ 

. 065822222776361  DO/ 
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172 

DATA 

Xt  102 

173 

DATA 

X(  103 

174 

DATA 

Xt  104 

175 

DATA 

Xt  105 

176 

DATA 

Xt  106 

177 

. 

DATA 

Xt  107 

178 

DATA 

Xt  108 

179 

DATA 

Xt  109 

180 

C 

N = 40 

181 

DATA 

Xt  110 

182 

DATA 

Xt  111 

183 

DATA 

Xt  112 

184 

DATA 

Xt  113 

185 

DATA 

Xt  114 

186 

DATA 

Xt  115 

187 

DATA 

Xt  116 

188 

DATA 

Xt  117 

189 

DATA 

Xt  118 

190 

DATA 

Xt  1.1-9 

191 

DATA 

Xt  120 

192 

DATA 

Xt  121 

193 

DATA 

Xt  122 

194 

DATA 

Xt  123 

195 

DATA 

Xt  124 

196 

DATA 

Xt  125 

197 

DATA 

Xt  126 

198 

DATA 

Xt  127 

199 

DATA 

Xt  128 

200 

DATA 

Xt  129 

201 

c 

N = 4B 

202 

DATA 

Xt  130 

203 

DATA 

X t 1 3 1 

204 

DATA 

Xt  132 

205 

DATA 

Xt  133 

206 

DATA 

Xt  134 

207 

DATA 

Xt  135 

208 

DATA 

Xt  136 

209 

DATA 

Xt  137 

210 

DATA 

Xt  138 

211 

DATA 

Xt  139 

212 

DATA 

Xt  140 

213 

DATA 

Xt  141 

214 

DATA 

Xt  142 

215 

DATA 

Xt  143 

216 

DATA 

Xt  144 

217 

DATA 

Xt  145 

218 

DATA 

Xt  146 

219 

DATA 

Xt  147 

220 

DATA 

Xt  148 

221 

DATA 

Xt  149 

222 

DATA 

Xt  150 

223 

DATA 

Xt  151 

224 

DATA 

Xt  152 

225 

DATA 

Xt  153 

226 

c 

N = 64 

227 

DATA 

Xt  154 

228 

DATA 

Xt  155 

)/0. 663044266930215D0/,  At 
)/0. 587715757240762D0/»  At 
)/0. 506899908932229D0//  At 
)/0. 4213512761  30635D0//  At 
)/0. 331868602282127D0/,  At 
) / O . 2392873622521 3700//  At 
>/0. 144471961582796D0/,  At 
)/0.048307665687738D0/,  At 

)/0. 998237709710559D0//  At 
)/0. 99072623B699457D0//  At 
) / O . 97  72  5 99  4 99  8 37  7’ 4 DO//  At 
)/0. 95791 68 19213791D0//  At 
)/0. 932812808278676D0//  At 
)/0. 902098806968874D0//  At 
)/0. 86595950321 225 9D0/»  At 
)/0. 82461223083331  IDO//  At 
)/0. 778305651426519D0//  At 
)/0. 727318255189927D0//  At 
)/0. 671956684614179D0/.  At 
)/0. 612553889667980D0/,  At 
)/0. 549467125095123D0//  At 
)/0. 483075801686178D0//  At 
)/0, 413779204371605D0//  At 
)/0. 341994090825758D0//  At 
)/0.268152185007253D0//  At 
)/0. 192697580701371D0//  At 
)/0. 1 16084O70675255D0//  At 
)/0. 03877241750605OD0//  At 

)/0. 998771007252426D0//  At 
)/0. 993530172266350D0/,  At 
)/0. 984124583722826D0/,  At 
) / O . 97059 1 59254 62 47D0//  At 
)/0. 952987703160430D0//  At 
)/0.931386690706554D0/,  At 
)/0. 905879136715569D0//  At 
)/0. e76572020274247D0//  At 
) /O . 84358826 162 43 93D0//  At 
)/0. 807066204029442D0/,  At 
)/0. 767159032515740D0//  At 
)/0.724034130923B14D0//  At 
)/0. 677872379632663D0//  At 
)/0. 628867396776513D0//  At 
)/0. 577224726083972D0//  At 
1/0.52316097 4 722233D0//  At 
1/0.466902904750958D0//  At 
)/0. 408686481990716D0//  At 
) / O . 34  8755886292 1 6 ODO//  At 
I/O. 2873624B7355455D0//  At 
1/0.224763790394689D0//  At 
I/O. 161222356068891D0//  At 
)/0. 097004699209462D0//  At 
)/0. 032380170962869D0//  At 

I/O. 999305041735772D0//  At 
1/0.9963401 1677 1955D0/,  At 


1021/0. 072345794 108B48D0/ 
1031/0.07819389578707000/ 
104 1/0. 08331 1924226946D0/ 
1051/0. 087652093004403D0/ 
1061/0. 091 173878695763D0/ 
1071/0. 093844399080804DO/ 
1081/0. 095638720079274D0/ 
1091/0. 0965400885 14727 DO/ 

1101/0. 00452 127709e533D0/ 
111  1/0. 010498284531 152D0/ 
1 121/0. 01 642105838 1907D0/ 
1131/0. 022245849 194 166D0/ 
1 14 1/0. 027937006980023D0/ 
1151/0. 033460 195282547D0/ 
1161/0. 038782 167974472D0/ 
1171/0. 0438709C8 1 85673D0/ 
1181/0. 048695807635072D0/ 
1191/0. 053227846983936D0/ 
1201/0. 05743976909939 IDO/ 
1211/0. 061 30624 2492928D0/ 
1221/0. 06480401 345660 IDO/ 
1231/0. 067912045815233D0/ 
124 1/0. 07061 1 64739 1286D0/ 
1251/0. 072886582395804D0/ 
1261/0. 074723169057968D0/ 
1271/0. 076 11 036 1900626D0/ 
1281/0. 07703981 8 164247 DO/ 
1291/0. 07 7 50594 7 97 8424 DO/ 

1301/0. 003 1 5334 6052305D0/ 
131 1/0. 007327553901 276D0/ 
1321/0. 01 1477234579234 DO/ 
1331/0. 01 55793 1572294 3D0/ 
1341/0. 0196161 60 4 57 35 5 DO/ 
1351/0. 023570760839324 DO/ 
1361/0. 027426509708356D0/ 
137 1 /O. 031 167227832798D0/ 
1381/0. 034777222564770D0/ 
1391/0. 03824 1351 065B30D0/ 
1401/0. 041545082943464D0/ 
141 1/0. 044674560856694D0/ 
1421/0. 04761 6658492490D0/ 
1431/0. 050359035553854 DO/ 
1441/0. 052890 189485 193D0/ 
1451/0. 055 199503699984D0/ 
1461/0. 057277292 100403D0/ 
1471/0. 059 11 4839698395D0/ 
1481/0. 060704439 165893D0/ 
1491/0. 062039423 159892D0/ 
1501/0. 0631 14 192286254 DO/ 
151 1/0. 063924238584648D0/ 
1521/0. 064466 164 435950D0/ 
1531/0. 06473769681 2683D0/ 

154 1/0. 00 178328072 1696D0/ 
1551/0. 004 14 7033260562D0/ 
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229 

DATA 

XC  156 

230 

DATA 

XC  157 

231 

DATA 

XC  158 

232 

DATA 

XC  159 

233 

DATA 

XC  160 

234 

DATA 

XC  161 

235 

DATA 

XC  162 

236 

DATA 

XC  163 

237 

DATA 

XC  164 

238 

DATA 

XC  165 

239 

DATA 

XC  166 

240 

DATA 

XC  167 

241 

DATA 

XC  168 

242 

DATA 

XC  169 

243 

DATA 

XC  170 

244 

DATA 

XC  171 

245 

DATA 

XC  172 

246 

DATA 

XC  173 

247 

DATA 

XC  174 

248 

DATA 

XC  175 

249 

DATA 

XC  176 

250 

DATA 

XC  177 

251 

DATA 

XC  178 

252 

DATA 

XC  179 

253 

DATA 

XC  180 

254 

DATA 

XC  181 

255 

DATA 

XC  182 

256 

DATA 

XC  183 

257 

DATA 

XC  184 

258 

DATA 

XC  185 

259 

C 

N = 80 

260 

DATA 

XC  186 

261 

DATA 

XC  187 

262 

DATA 

XC  188 

263 

DATA 

XC  189 

264 

DATA 

XC  190 

265 

DATA 

XC  191 

266 

DATA 

XC  192 

267 

DATA 

XC  193 

268 

DATA 

XC  194 

269 

DATA 

XC  195 

270 

DATA 

XC  196 

271 

DATA 

XC  197 

272 

DATA 

XC  198 

273 

DATA 

XC  199 

274 

DATA 

X C 200 

275 

DATA 

X C 201 

276 

DATA 

X C 202 

277 

DATA 

X C 203 

278 

DATA 

X C 204 

279 

DATA 

X C 205 

280 

DATA 

X C 206 

281 

DATA 

XC  207 

282 

DATA 

X ( 208 

283 

DATA 

X C 209 

284 

DATA 

XC210 

285 

DATA 

X C 21  1 

)/0. 991013371476744D0/ 
)/0. 983336253884625D0/ 
) / 0 . 9 7 33268 2 77 8 991 ODO/ 
>/0. 961008799652053D0/ 
)/0. 94641137  4858402D0/ 
)/0. 929569 1721 3 1939D0/ 
)/0. 9 105221 37078502D0/ 
)/0. 8893154459951 14D0/ 
I/O. 865999398 154092D0/ 
)/0. 840629296252580D0/ 
)/0 . 8132653 15122797D0/ 
)/0. 78 39723589 4 334 IDO/ 
)/0. 7528 1990726053 IDO/ 
)/0. 71 98818501 7 1610D0/ 
)/0. 68523631 3054233D0/ 
)/0 . 64896547 1254657D0/ 
)/0.611155355172393D0/ 
)/0.571B95646202634DO/ 
)/0. 53 12794640 19894D0/ 
)/0. 4894031 45707052D0/ 
)/0. 446366017253464D0/ 
)/0. 4022701 5796399 IDO/ 
)/0. 357220 158337668D0/ 
)/0. 3 1132287 199021 ODO/ 
)/0,264687162208767D0/ 
I/O. 2 1742364 3740007D0/ 
)/0. 169644420423992D0/ 
)/0. 1214628 19296 120D0/ 
)/0. 0729931 2 1787799D0/ 
)/0. 024350292663424D0/ 

)/0. 99955382265 1630D0/ 
)/0 . 99764S864398237D0/ 
)/0.994227540S65688D0/ 
)/0. 98929 1302499755D0/ 
)/0.982848572738629D0/ 
)/0. 974909 140585727D0/ 
)/0. 965485089043799D0/ 
)/0. 95459076634363400/ 
)/0.942242761309872D0/ 
)/0. 928459877 172445D0/ 
)/0 . 9132631 02571757D0/ 
)/0. 896675579438770D0/ 
)/0. 87872256767821 3D0/ 
)/0. 85943 140666311 IDO/ 
)/0. 83883 147 3580255D0/ 
)/0. 8 169541 3868 1463D0/ 
)/0. 7938327 17504605D0/ 
)/0. 7695024  20 13504  IDO/ 
)/0. 744000297583597D0/ 
)/0. 7173651B5362099D0/ 
)/0. 68963764434 2027D0/ 
)/0. 6608598989861 19D0/ 
)/0. 63 107577304 687 IDO/ 
)/0. 60033062282975 IDO/ 
)/0. 56867 1 2681 22709D0/ 
)/0. 5361 4 59208971 31  DO/ 


AC 156)/0. 
A ( 1 57 ) /O . 
A ( 1 58 ) /O . 
A < 1 59 ) /O . 
AC  1601/0. 
AC  161 1/0. 
AC  1621/0. 
AC  1631/0. 
AC  1641/0. 
AC  1651/0. 
AC  1661/0. 
AC  1671/0. 
AC  1681/0. 
AC  1691/0. 
AC  1701/0. 
AC  171 1/0. 
AC  1721/0. 
AC  1731/0. 
AC1741/0. 
AC  1751/0. 
AC  1761/0. 
AC  1771/0. 
AC  1781/0. 


A 

C 179 

)/0. 

A 

C 180 

)/0. 

A 

C 181 

O 

\ 

/A 

A 

C 182 

)/0. 

A 

C 183 

)/0. 

A 

C 184 

I/O. 

A 

C 185 

)/0. 

A 

C 186 

)/0. 

A 

C 187 

)/0. 

A 

c ie8 

)/0. 

A 

C 189 

)/0. 

A 

C 190 

)/0. 

A 

C 191 

)/0. 

A 

C 192 

)/0. 

A 

C 193 

>/0. 

A 

C 194 

)/0. 

A 

C 195 

)/0. 

A 

C 196 

)/0. 

A 

C 197 

)/0. 

A 

C 198 

)/0. 

A 

C 199 

)/0. 

A 

C 200 

)/0. 

A 

C 201 

)/0. 

A 

C 202 

)/0. 

A 

C 203 

)/  0. 

A 

C 204 

>/0. 

A 

C 205 

)/0. 

A 

C 206 

)/0. 

A 

C 207 

)/0. 

A 

C 208 

)/0. 

A 

C 209 

)/0. 

A 

C 210 

)/0. 

A 

C 2 1 1 

)/0. 

006504 457 96897 8D0/ 
00884 6759826363D0/ 
011168139460131D0/ 
01 34630478967 18D0/ 
01 572603047 6024D0/ 
0179517 1 5 77 5697 DO/ 
020 1 34823 153530D0/ 
022270 173808383D0/ 
0243527025687 10D0/ 
026 37 7 4697 15 054 DO/ 
02833967261 4259D0/ 
030234657072402D0/ 
03205792835485 IDO/ 
033805161 837 141D0/ 
035472213256882D0/ 
037055 128540240D0/ 
038550 1 53 1 7861 5D0/ 
03995374 1132720D0/ 
041 26 2 56 324 2 62 3 DO/ 
042473515123653D0/ 
04 358 372 4 52 9 32 3 DO/ 
044590558 1 63756D0/ 
045491 6279274 18D0/ 
046284796581 31 4D0/ 
046968 1828 1 621 ODO/ 
047540 1 657 14B30D0/ 
047999388596458D0/ 
04 8 34 4 76 2 23 4 80 2 DO/ 
04857546744 1503D0/ 
048690957009 139D0/ 

001 1 44950003 186D0/ 
00266353358951 2D0/ 
004 180313124694D0/ 
00569092245 1403D0/ 
0071 92904768 117D0/ 
00 8683 94 526 9 26 ODO/ 
010161766041 103D0/ 
011624114120797D0/ 
01 306876 159240 IDO/ 
014493508040509D0/ 
01 5896 183583725 DO/ 
017274652056269D0/ 
01 862681 4208299D0/ 
019950610878141D0/ 
021 244026 115782D0/ 
022505090246332D0/ 
02373 1B82865930D0/ 
0249225357641 15D0/ 
026075235767565D0/ 
027188227500486DO/ 
02825981 6057276D0/ 
029288369583267D0/ 
03027232 1 759557D0/ 
031210174188114D0/ 
0321 00498673487 DO/ 
03294 1939397645D0/ 


286 

DATA 

X ( 21 2 

287 

DATA 

X ( 21 3 

288 

DATA 

X ( 214 

289 

DATA 

X ( 215 

290 

DATA 

X ( 216 

291 

DATA 

X ( 217 

292 

DATA 

X(  218 

293 

DATA 

X ( 219 

294 

DATA 

X ( 220 

295 

DATA 

X ( 221 

296 

DATA 

X ( 222 

297 

DATA 

X ( 223 

298 

DATA 

X ( 224 

299 

DATA 

X ( 225 

300 

C 

N = 96 

301 

DATA 

X ( 226 

302 

DATA 

X ( 227 

303 

DATA 

X ( 228 

304 

DATA 

X ( 229 

305 

DATA 

X ( 230 

306 

DATA 

X ( 231 

307 

DATA 

X ( 232 

308 

DATA 

X ( 233 

309 

DATA 

X ( 234 

310 

DATA 

X ( 235 

311 

DATA 

X ( 236 

312 

DATA 

X ( 237 

313 

DATA 

X ( 238 

314 

DATA 

X ( 239 

315 

DATA 

X ( 240 

316 

DATA 

X ( 24 1 

317 

DATA 

X ( 242 

318 

DATA 

X ( 243 

319 

DATA 

X ( 244 

320 

DATA 

X ( 245 

321 

DATA 

X ( 246 

322 

DATA 

X ( 247 

323 

DATA 

X ( 248 

324 

DATA 

X ( 249 

325 

DATA 

X ( 250 

326 

DATA 

X ( 251 

327 

DATA 

X ( 252 

328 

DATA 

X ( 253 

329 

DATA 

X ( 254 

330 

DATA 

X ( 255 

331 

DATA 

X ( 256 

332 

DATA 

X ( 257 

333 

DATA 

X ( 258 

334 

DATA 

X ( 259 

335 

DATA 

X ( 260 

336 

DATA 

X ( 261 

337 

DATA 

X ( 262 

338 

DATA 

X ( 263 

339 

DATA 

X ( 264 

340 

DATA 

X ( 265 

341 

DATA 

X ( 266 

342 

DATA 

X ( 267 

)/0. 5028041 118887B4D0/ 
)/0. 4686966 151 70544D0/ 
)/0 . 4 3387537 0831 7 56D0/ 
>/0.39B393405881969DO/ 
>/0. 362304753499487D0/ 
)/0. 32566437074770100/ 
)/0. 28852805488451 100/ 
) /0 . 25 095235 839 227 2D0/ 
) /0 . 212994502857666D0/ 
)/0. 174712291832646D0/ 
)/0. 136164022809143D0/ 
)/0. 09740839844 1584D0/ 
)/0. 058504437152420DO/ 
)/0. 01951 13B3256793D0/ 

)/0 . 99 968950 38 83230D0/ 
)/0. 9983643758631 81 DO/ 
)/0. 99598 1842987209D0/ 
) /0 . 9 92 5 4 390 03 2 37 62D0/ 
)/0. 988054 126329623D0/ 
)/0 . 9 8 25 1726 356 30 14 DO/ 
)/0. 9759391 745851 36D0/ 
)/0. 96832682846326400/ 
)/0. 95968829 1448742D0/ 
)/0. 950032717784437D0/ 
)/0 . 9 3937 03397 5 27 5 5D0/ 
)/0. 9277 124567 223 08D0/ 
)/0. 9 1507 14231 20898D0/ 
)/0. 901460635315852D0/ 
)/0. 88689451740242000/ 
)/0. 87138850590929600/ 
)/0. 8549590334 3460 IDO/ 
)/0. 83762351 12281 87D0/ 
)/0. 81 94003 1073793 IDO/ 
)/0. 8003087441 391 40D0/ 
) / O . 7 8036904 38674 3 3D0/ 
)/0. 75960234 11 76647D0/ 
) / O . 7 3803064 374 44 OODO/ 
)/0. 7 1567681 234 8967D0/ 
)/0. 69256453664217100/ 
)/0. 6687 183 1 00439 16D0/ 
) /O . 644 1 634 037 8496 7D0/ 
)/0. 6 189258401 254 68D0/ 
)/0. 593032364777572D0/ 
)/0. 5665 104 18561 397D0/ 
)/0. 5393881 08324357D0/ 
) /O. 51 16941 771 54667D0/ 
)/0. 483457973920596D0/ 
)/0. 454709422167743D0/ 
)/0. 425478988407300D0/ 
)/0. 395797649828908D0/ 
)/0. 365696861472313D0/ 
)/0. 335208522892625D0/ 
)/0. 304364S4435449600/ 
)/0.273198B12591049D0/ 
)/0. 24 17431 561 63840D0/ 
)/0. 21 00313 10460567D0/ 


A ( 212 )/0. 03373321490461  IDO/ 
A<213)/0. 034473 120451 75 3D0/ 
A( 214 ) /O . 035 16 0529 04 4 74 7 DO/ 
A ( 215 )/0. 0357943939534 16D0/ 
A(216)/0. 03637374  990 5835 DO/ 
A(217)/0. 0368977 14638276D0/ 
A(218)/0. 037365490238730D0/ 
A<219)/0. 037776364  362001  DO/ 
A < 220 )/0. 038 12971 1314 477D0/ 
A(221 )/0. 038424 99 300 695 9 DO/ 
A(222)/0. 038 66 175 977 4 07 6 DO/ 
A ( 223 )/0. 03883965 105905 IDO/ 
A<224)/0. 038958395962769D0/ 
A ( 225 ) /O . 03901781 3656306D0/ 

A ( 226 ) /O . 000796792065552D0/ 
A(227 )/0. 00 185 3960 78894 6D0/ 
A( 228 )/0. 00291 073 181 7934D0/ 
A( 229 )/0. 003964554338444 DO/ 
A ( 230 ) /O . 00501 4 20 274 2 927 DO/ 
A( 231 )/0. 006058545504235D0/ 
A( 232 )/0. 007096470791 153D0/ 
A ( 233 ) /O . 00812687  6925698D0/ 
A<234)/0. 009 1 4067 1230783D0/ 
A ( 235 ) /O .010 1 6 077 053 5008 DO/ 
A ( 236 ) /O . 01 1 162102099838D0/ 
A(237)/0. 012151604671088D0/ 
A ( 238 )/0. 01 3 12822956696 IDO/ 
A < 239 )/0. 01 409094 17723 14D0/ 
A<  240) /O. 01 503872 1026994D0/ 
A ( 241 ) /O . 015970562902562D0/ 
A ( 242 )/0. 01 688547986424 5D0/ 
A ( 243 )/0. 01 77825023 16045D0/ 
A ( 244 )/0. 01866067962741  IDO/ 
A ( 245 )/0. 019519081 1401 45D0/ 
A ( 246 )/0. 020356797 154 333D0/ 
A(247)/0. 02 1172S39892191D0/ 
A(248)/0. 02 196664 4438744D0/ 
A ( 249 ) /O . 022737069658329D0/ 
A ( 250 ) /O . 0234 8 3 39 908 5 926 DO/ 
A ( 25 1 )/0. 02420484 1792364 DO/ 
A(252)/0. 02 4 90063 3 22248 3D0/ 
A( 253 )/0. 02557003600534 9D0/ 
A( 254 )/0. 02621 2340735672D0/ 
A C 255 )/0. 02682686672559  IDO/ 
A C 256 ) /O . 027 412962726029D0/ 
A ( 257 ) /G . 02797000761 684 8D0/ 
A ( 258 )/0. 02849741 1065085D0/ 
A ( 259 ) /O . 02899461 4 150555D0/ 
A( 260) /O. 02946 1 089958 167D0/ 
A (261 )/0. 029896344 136328D0/ 
A(262)/0. 0302999 15420827D0/ 
A ( 263 )/0. 03067 1376 123669D0/ 
A ( 264 )/0. 03 101 033258631 3D0/ 
A ( 265 ) /0 . 031316425596861  DO/ 
A(266)/0. 031 589330 77 0727 DO/ 
A ( 267 ) /O . 03182875889441  IDO/ 


85 


343 

DATA  X(26B)/0. 17809688236761800/# 

A ( 268  )/0. 032034 45623 1992D0/ 

344 

DATA  XC269)/0. 1 459737 1 4654896D0/# 

AC  269 )/0. 032206204  794  030D0/ 

345 

DATA  XC270)/0. 1 1 36958501 1 0665D0/# 

AC  270 )/0. 03234 382256B575D0/ 

346 

DATA  XC271 )/0.081297495464425D0/# 

AC271 )/0. 032447 1 637 14064D0/ 

347 

DATA  X<272)/0. 048812985136049D0/# 

AC272)/0. 03251611871 3 86 8 DO/ 

348 

DATA  X(273)/0. 016276744849602D0/# 

AC273)/0. 03255061 4492363D0/ 

349 

C 

350 

C 

Test  N 

351 

N = NX 

352 

ALPHA=0 . 5 DO* C BX+AX ) 

353 

BETArO. 5D0*CBX-AX) 

354 

IF  CN.EQ.l)  THEN 

355 

Z( 1 ) = ALPHA 

356 

W(  1 ) rBX-AX 

357 

ELSE  IF  ( (N.GT. 1 . AND. N. LE . 16) .OR. N 

. EQ.20.0R.N.EQ.24.0R.N.EQ.32.0R. 

358 

+ N. EQ. 40. OR. N. EQ. 48. OR. N. EQ. 64 . OR 

.N.EQ.80.0R.N.EQ.96)  THEN 

359 

c 

Set  K equal  to  initial  subscript  and  store  results 

360 

KrKTABCN) 

361 

M=N/2 

362 

DO  904  J=1,M 

363 

JTAB=K-1+J 

364 

WTEMP=BETA*A( JTAB) 

365 

DELTArBETA*X( JTAB) 

366 

Z(J)= ALPHA -DELTA 

367 

W( J) =WTEMP 

368 

JPrN+l-J 

369 

Z( JP)=ALPHA+DELTA 

370 

W( JP)=WTEMP 

371 

904 

CONTINUE 

372 

IF  ( (N-M-M) . EQ. 0)  RETURN 

373 

Z ( M+ 1 ) = ALFHA 

374 

JMID=K+M 

375 

UKM  + l 5 =BETA*A(  JMID) 

376 

ELSE 

377 

ZN  = N 

378 

WRITE  ( 1 » 50 ) ZN 

379 

50 

FORMAT( 'OError  in  GSET . N has 

the  nonperm i ss ab 1 e value'# 

380 

+ 1 PEI 1.3) 

381 

ENDIF 

382 

RETURN 

383 

END 
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1 

STITL 

EUALI  --  61-point  Gauss-Kronrod  integration 

2 

FUNCTION  EUALI (A, B, FUNC, NS) 

3 

C 

Function  to  evaluate  the  integral  of  the  function  FUNC 

4 

C 

between  A and  B.  The  number  of  intervals  starts  at  NSTEP 

5 

C 

as  defined  -in  EUICOM  and  then  goes  up  by  a factor  of  two 

6 

C 

until  the  error  in  the  integral  is  smaller  than  EPS. 

7 

C 

In  this  version/  61-point  Gauss-Kronrod  integration 

8 

C 

using  the  subroutine  DQK61  from  the  NBS  library  CMLIB 

9 

C 

is  done. 

10 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

11 

C 

The  maximum  number  of  times  that  the  integral  is  evaluated 

12 

PARAMETER  (IDIUMX=10) 

13 

EXTERNAL  FUNC 

14 

COMMON  /EUICOM/  NSTEP/ EPS 

15 

NSrNSTEP 

16 

IDIUrO 

17 

c 

18 

c 

Start  of  loop 

19  * 

4 

STEP: (B-A)/NS 

20 

EU=0. DO 

21 

ERR;0. DO 

22 

DO  902  1:1/ NS 

23 

AX : A+ ( 1-1 ) *STEP 

24 

BX : AX+STEP 

25 

CALL  DQK61 (FUNC/ AX/ BX/ RESULT/ ABSERR / RESABS / RESASC ) 

26 

EU : EU  + RESUL  T 

27 

ERR:ERR+ABSERR 

28 

902 

CONTINUE 

29 

IF  (EPS.LE.O. DO)  THEN 

30 

CD 

WRITE  (1,*)  A, B, ERR, DABS(EU*EPS) 

31 

IF  (ERR. LE. DABS(EU*EPS) ) GO  TO  6 

32 

ELSE 

33 

CD 

WRITE  (1/*)  A, B, ERR, EPS 

34 

IF  (ERR.LE.EPS)  GO  TO  6 

35 

ENDIF 

36 

IF  (NS. NE. NSTEP)  THEN 

37 

IF  (EPS. LE. 0. DO)  THEN 

38 

IF  ( DABS ( EU-EUAL 10) . LE . DABS ( EU* EPS ) ) GO  TO  6 

39 

CD 

WRITE  (1,*)  A , B, DABS ( EU-EUAL 10 ) , DABS ( EU*EPS ) 

40 

ELSE 

41 

IF  (DABS(EU-EUALIO) . LE. EPS)  GO  TO  6 

42 

CD 

WRITE  (1/*)  A, B, DABS (EU-EUALIO) , EPS 

43 

ENDIF 

44 

ENDIF 

45 

EUAL 10; EU 

46 

NS:2*NS 

47 

I DI U = I D IU+ 1 

48 

IF  ( IDIU. LE. IDIUMX)  GO  TO  4 

49 

NS:  — 1 

50 

6 

EUAL I : EU 

51 

RETURN 

52 

END 
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1 

2 

3 

4 

5 

6 

7 

Q 

9 

10 

1 1 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 


DIU.  - K.U.LEUUEN 
DIO.  - K . U . LEUUEN 


SUBROUTINE  DQK6 1 C F,  A, B, RESULT, ABSERR, RESABS , RES A SC ) 
C***BEG I N PROLOGUE  DQK61 
C***DATE  WRITTEN  800101  (YYMMDD) 

C***REUISION  DATE  810101  (YYMMDD) 

C*** CATEGORY  NO.  D1 

C* <K<<KEYWORDS  61-POINT  GAUSS-KRONROD  RULES 
C***AUTHOR  PIESSENS, ROBERT, APPL.  MATH.  & PROGR. 

C DE  DONCKER, ELISE,APPL.  MATH.  & PROGR 

C D.K.  KAHANER  ( NBS ) 

C**#PURPOSE  TO  COMPUTE  I r INTEGRAL  OF  F OUER  <A,B)  WITH  ERROR 
C ESTIMATE 

C J = INTEGRAL  OF  DABS ( F ) OUER  ( A , B ) 

c***description 

c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


INTEGRATION  rule 

STANDARD  FORTRAN  SUBROUTINE 

DOUDLE  PRECISION  UERSION 


PARAMETERS 
ON  ENTRY 
F 


DOUBLE  PRECISION 

SUBROUTINE  DEFINING  THE  INTEGRAND  F(X)  AT  X: 

CALL  F(X,F(X>).  THE  ACTUAL  NAME  FOR  F NEEDS  TO  BE 
DECLARED  EXTERNAL  IN  THE  CALLING  PROGRAM. 

DOUBLE  PRECISION 

LOWER  LIMIT  OF  INTEGRATION 

DOUBLE  PRECISION  • 

UPPER  LIMIT  OF  INTEGRATION 


DOUBLE  PRECISION 

APPROXIMATION  TO  THE  INTEGRAL  I 
RESULT  IS  COMPUTED  BY  APPLYING  THE  61-POINT 
KRONROD  RULE  (RESK)  OBTAINED  BY  OPTIMAL  ADDITION  OF 
ABSCISSAE  TO  THE  30-POINT  GAUSS  RULE  CRESG). 


ABSERR  - DOUBLE  PRECISION 

ESTIMATE  OF  THE  MODULUS  OF  THE  ABSOLUTE  ERROR, 
WHICH  SHOULD  EQUAL  OR  EXCEED  DABS ( I -RESULT ) 

RESABS  - DOUBLE  PRECISION 

APPROXIMATION  TO  THE  INTEGRAL  J 


ON  RETURN 
RESULT  - 


RESASC  - 


DOUBLE  PRECISION 

APPROXIMATION  TO  THE  INTEGRAL  OF  DABS ( F- 1/ ( B-A ) ) 


♦♦♦REFERENCES  (NONE) 

♦ ♦♦ROUTINES  CALLED  D 1 MACH , DABS, DMAX 1 , DMI N 1 , F ( USER-PROUI DED  FUNCTION) 


♦♦♦END  PROLOGUE  DQK61 
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58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 


DOUBLE  PRECISION  A, DABSC, ABSERR, B,  CENTR, DABS, DHLGTH, DMAX1 , DM  INI , 

* D1MACH, EPMACH, FC, FSUM, FUAL1,  FUAL2, FU1 , FU2, HLGTH, RESABS, RESASC- 

* RESG, RESK, RESKH, RESULT , UFLOW, WG, WGK, XGK 
INTEGER  J , JTW » JTWM 1 

DIMENSION  FU1 ( 30 ) , FU2 ( 30 ) , XGK< 31 ) , WGK ( 31 ) , WG C 1 5 ) 

THE  ABSCISSAE  AND  WEIGHTS  ARE  GIUEN  FOR  THE 

INTERUAL  (-1,1).  BECAUSE  OF  SYMMETRY  ONLY  THE  POSITIVE 

ABSCISSAE  AND  THEIR  CORRESPONDING  WEIGHTS  ARE  GIUEN. 

XGK  - ABSCISSAE  OF  THE  61-POINT  KRONROD  RULE 

XGK ( 2 ) , XGK ( 4 ) ...  ABSCISSAE  OF  THE  30-POINT 

GAUSS  RULE 

XGK ( 1 ) , XGK ( 3 ) ...  OPTIMALLY  ADDED  ABSCISSAE 

TO  THE  30-POINT  GAUSS  RULE 

WGK  - WEIGHTS  OF  THE  61-POINT  KRONROD  RULE 

WG  - WEIGTHS  OF  THE  30-POINT  GAUSS  RULE 


GAUSS  QUADRATURE  WEIGHTS  AND  KRONRON  QUADRATURE  ABSCISSAE  AND  WEIGHTS 


81 

C 

AS  EVALUATED 

WITH 

80 

DECIMAL  DIGIT  ARITHMETI 

C BY  L.  W. 

FULLERTON 

82 

C 

BELL  LABS,  NOV. 

1981  . 

83 

C 

84 

DATA 

WG 

( 

1 ) 

✓ 

0.0079681924 

9616660561 

5465883474 

674 

DO 

/ 

85 

DATA 

WG 

( 

2) 

/ 

0.0184664683 

1109095914 

2302131912 

047 

DO 

/ 

86 

DATA 

WG 

( 

3) 

s 

0.0287847078 

8332336934 

9719179611 

292 

DO 

/ 

87 

DATA 

WG 

( 

4) 

/ 

0.0387991925 

6962704959 

6801936446 

348 

DO 

/ 

88 

DATA 

WG 

( 

5) 

s 

0.0484026728 

3059405290 

2938140422 

808 

DO 

/ 

89 

DATA 

WG 

( 

6) 

/ 

0. 057.4931562 

1761906648 

1721689402 

056 

DO 

/ 

90 

DATA 

WG 

< 

7) 

/ 

0.0659742298 

8218049512 

8128515115 

962 

DO 

/ 

91 

DATA 

WG 

( 

8) 

/ 

0.0737559747 

3770520626 

8243850022 

191 

DO 

/ 

92 

DATA 

WG 

( 

9) 

/ 

0.0807558952 

2942021535 

4694938460 

530 

DO 

/ 

93 

DATA 

WG 

( 

10) 

s 

0.0868997872 

0108297980 

2387530715 

126 

DO 

/ 

94 

DATA 

WG 

( 

11  ) 

/ 

0. 0921225222 

3778612871 

7632707087 

619 

DO 

/ 

95 

DATA 

WG 

( 

12) 

/ 

0.0963687371 

7464425963 

9468626351 

810 

DO 

/ 

96 

DATA 

WG 

( 

13) 

s 

0.0995934205 

8679526706 

2780282103 

569 

DO 

/ 

97 

DATA 

WG 

( 

14) 

s 

0. 1017623897 

4840550459 

6428952168 

554 

DO 

✓ 

98 

DATA 

WG 

( 

15) 

/ 

0. 1028526528 

9355884034 

1285636705 

415 

DO 

/ 

99 

C 

100 

DATA 

XGK 

( 

1 ) 

/ 

0.9994844100 

5049063757 

1 3258S5705 

811 

DO 

/ 

101 

DATA 

XGK 

( 

2) 

/ 

0.9968934840 

7464954027 

1630050918 

695 

DO 

/ 

102 

DATA 

XGK 

< 

3) 

/ 

0.9916309968 

7040459485 

8628366109 

486 

DO 

/ 

103 

DATA 

XGK 

< 

4) 

s 

0. 9836681232 

7974720997 

0032581605 

663 

DO 

/ 

104 

DATA 

XGK 

< 

5) 

/ 

0.9731163225 

0112626837 

4693868423 

707 

DO 

✓ 

105 

DATA 

XGK 

( 

6) 

/ 

0.9600218649 

6830751221 

6871025581 

798 

DO 

✓ 

106 

DATA 

XGK 

( 

7) 

/ 

0. 9443744447 

4855997941 

5831324037 

439 

DO 

/ 

107 

DATA 

XGK 

( 

8) 

s 

0.9262000474 

2927432587 

9324277080 

474 

DO 

/ 

108 

DATA 

XGK 

( 

9) 

/ 

0.9055733076 

9990779854 

6522558925 

958 

DO 

/ 

109 

DATA 

XGK 

( 

10) 

/ 

0. 8825605357 

9205268154 

31 16462530 

226 

DO 

/ 

110 

DATA 

XGK 

( 

11  ) 

/ 

0.8572052335 

4606109895 

8658510658 

944 

DO 

/ 

111 

DATA 

XGK 

( 

12) 

/ 

0.8295657623 

8276839744 

28981 19732 

502 

DO 

/ 

112 

DATA 

XGK 

< 

13) 

s 

0. 7997278358 

2183908301 

3668942322 

683 

DO 

/ 

1 1 3 

DATA 

XGK 

( 

14  ) 

/ 

0. 7677774321 

0482619491 

797734097 4 

503 

DO 

/ 

114 

DATA 

XGK 

( 

15) 

/ 

0.7337900624 

5322680472 

6171131369 

528 

DO 

/ 

89 


115 

DATA 

XGK 

( 

16) 

/ 

0.6970504947 

9331579693 

2292388026 

640 

DO 

/ 

116 

DATA 

XGK 

( 

17) 

/ 

0. 6600610641 

2662696137 

0053668149 

271 

DO 

/ 

117 

DATA 

XGK 

< 

18) 

✓ 

0.6205261029 

8924286114 

0477556431 

109 

DO 

/ 

118 

DATA 

XGK 

( 

19) 

/ 

0.5793452358 

2636169175 

6024932172 

540 

DO 

/ 

119 

DATA 

XGK 

( 

20) 

/ 

0.5366241481 

4201989926 

4169793311 

073 

DO 

/ 

120 

DATA 

XGK 

( 

21 ) 

/ 

0.4924804678 

6177857499 

3693061207 

709 

DO 

/ 

121 

DATA 

XGK 

( 

22) 

/ 

0.4470337695 

3808917678 

0609900322 

054 

DO 

/ 

122 

DATA 

XGK 

< 

23) 

/ 

0.4004012548 

3039439253 

5476211542 

661 

DO 

/ 

123 

DATA 

XGK 

< 

24) 

/ 

0.3527047255 

3087011347 

1037207089 

374 

DO 

/ 

124 

DATA 

XGK 

( 

25) 

/ 

0.3040732022 

7362507737 

2677107199 

257 

DO 

/ 

125 

DATA 

XGK 

( 

26) 

/ 

0.2546369261 

6788984643 

9805129817 

805 

DO 

/ 

126 

DATA 

XGK 

( 

27) 

/ 

0.2045251166 

8230909143 

8957671002 

025 

DO 

/ 

127 

DATA 

XGK 

( 

28) 

/ 

0. 1530699136 

0858354696 

3794672743 

256 

DO 

/ 

120 

DATA 

XGK 

( 

29) 

/ 

0. 1028069379 

6673703014 

7096751318 

001 

DO 

/ 

129 

DATA 

XGK 

( 

30) 

/ 

0.0514718425 

5531769583 

3025213166 

723 

DO 

/ 

130 

DATA 

XGK 

( 

31  ) 

/ 

0.0000000000 

0000000000 

0000000000 

000 

DO 

/ 

131 

C 

132 

DATA 

WGK 

( 

1) 

/ 

0.0013890136 

9867700762 

4551591226 

760 

DO 

/ 

133 

DATA 

WGK 

( 

2) 

/ 

0.0038904611 

2709988405 

1267201844 

516 

DO 

/ 

134 

DATA 

WGK 

< 

3) 

/ 

0.0066307039 

1593129217 

3319826369 

750 

DO 

/ 

135 

DATA 

WGK 

( 

4) 

/ 

0.0092732796 

5951776342 

8441146892 

024 

DO 

/ 

136 

DATA 

WGK 

( 

5) 

/ 

0.0118230152 

5349634174 

2232898853 

251 

DO 

/ 

137 

DATA 

WGK 

< 

6) 

/ 

0.0143697295 

0704580481 

2451432443 

580 

DO 

/ 

138 

DATA 

WGK 

( 

7) 

/ 

0.0169208891 

8905327262 

7572289420 

322 

DO 

/ 

139 

DATA 

WGK 

( 

8) 

/ 

0.0194141411 

9394238117 

3408951050 

128 

DO 

/ 

140 

DATA 

WGK 

< 

9) 

/ 

0.0210280358 

2160919229 

7167485738 

339 

DO 

/ 

141 

DATA 

WGK 

( 

10) 

/ 

0.0241911620 

7808060136 

5686370725 

232 

DO 

/ 

142 

DATA 

WGK 

( 

11  ) 

/ 

0.0265099540 

8233310161 

0601709335 

075 

DO 

/ 

143 

DATA 

WGK 

( 

12) 

/ 

0.0287540487 

6504129284 

3978785354 

334 

DO 

/ 

144 

DATA 

WGK 

( 

13) 

/ 

0.0309072575 

6238776247 

2884252943 

092 

DO 

/ 

145 

DATA 

WGK 

< 

14) 

/ 

0.0329814470 

5740372603 

1814191016 

854 

DO 

/ 

146 

DATA 

WGK 

( 

15) 

/ 

0.0349793380 

2806002413 

7499670731 

468 

DO 

/ 

147 

DATA 

WGK 

( 

16) 

/ 

0.0368823646 

5182122922 

3911065617 

136 

DO 

/ 

148 

DATA 

WGK 

C 

17) 

/ 

0.0386789456 

2472759295 

0348651532 

281 

DO 

/ 

149 

DATA 

WGK 

< 

18) 

/ 

0.0403745389 

515359591 1 

1995279752 

468 

DO 

/ 

150 

DATA 

WGK 

( 

19) 

/ 

0.0419698102 

1516424614 

7147541205 

970 

DO 

/ 

151 

DATA 

WGK 

( 

20) 

/ 

0. 0434525397 

0135606931 

6831728117 

073 

DO 

/ 

152 

DATA 

WGK 

< 

21  ) 

/ 

0.0448148001 

3316266319 

2355551616 

723 

DO 

/ 

153 

DATA 

WGK 

( 

22) 

/ 

0.0460592382 

7100698811 

6271735559 

374 

DO 

/ 

154 

DATA 

WGK 

( 

23) 

/ 

0.0471855465 

6929915394 

5261470181 

099 

DO 

/ 

155 

DATA 

WGK 

( 

24) 

/ 

0.0401858617 

5708712914 

0779492298 

305 

DO 

/ 

156 

DATA 

WGK 

( 

25) 

/ 

0.0490554345 

5502977888 

7528165367 

238 

DO 

/ 

157 

DATA 

WGK 

c 

26) 

/ 

0.0497956834 

2707420635 

7811569379 

942 

DO 

/ 

158 

DATA 

WGK 

( 

27) 

/ 

0. 0504059214 

0278234684 

0893085653 

585 

DO 

/ 

159 

DATA 

WGK 

< 

28) 

/ 

0.0508817958 

9874960649 

2297473049 

805 

DO 

/ 

160 

DATA 

WGK 

( 

29) 

/ 

0.0512215478 

4925877217 

0656282604 

944 

DO 

/ 

161 

DATA 

WGK 

( 

30) 

/ 

0.0514261285 

3745902593 

3862879215 

781 

DO 

/ 

162 

DATA 

WGK 

( 

31  ) 

/ 

0.0514947294 

2945156755 

8340433647 

099 

DO 

/ 

163 

C 

164 

C 

LIST 

OF 

MAJOR  VARIABLES 

165 

C 

166 

C 

167 

C 

CENTR 

MID  POINT  OF  THE  INTERVAL 

168 

C 

HLGTH 

1 - 

HALF-LENGTH  OF 

THE  INTERVAL 

169 

c 

DABSC 

ABSCISSA 

170 

c 

FVAL* 

FUNCTION  VALUE 

171 

c 

RESG 

- 

RESULT  OF  THE  30-POINT  GAUSS  RULE 

90 


172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 


RESK  - RESULT  OF  THE  61-POINT  KRONROD  RULE 
RESKH  - APPROXIMATION  TO  THE  MEAN  VALUE  OF  F 
OVER  < A , B ) , I.E.  TO  I/(B-A) 

MACHINE  DEPENDENT  CONSTANTS 


EPMACH  IS  THE  LARGEST  RELATIUE  SPACING. 
UFLOW  IS  THE  SMALLEST  POSITIVE  MAGNITUDE. 


EPMACH  = D1MACH ( 4 ) 

UFLOW  = DlMACH(l) 

C 

CENTR  = O. 5D+00*<B+A) 

HLGTH  = 0. 5D+00*(B-A) 

DHLGTH  = DABS ( HLGTH ) 

C 

C COMPUTE  THE  61-POINT  KRONROD  APPROXIMATION  TO  THE 

C INTEGRAL,  AND  ESTIMATE  THE  ABSOLUTE  ERROR. 

C 

C***FIRST  EXECUTABLE  STATEMENT  DQK61 


193 

RESG  = O.OD+OO 

194 

CALL  F < CENTR, FC  ) 

195 

RESK  r WGK(31)*FC 

1 196 

RESABS  = DABS ( RESK ) 

197 

DO  10  J=l, 15 

198 

JTW  = J*2 

199 

DABSC  = HLGTH*XGK( JTW) 

200 

CALL  F ( CENTR-DABSC, FUAL1 ) 

201 

CALL  F ( CENTR+ DABSC, FOAL 2 ) 

202 

FOl(JTW)  = F0AL1 

203 

F02 (JTW ) = F0AL2 

204 

FSUM  = FOAL 1 +F0AL2 

205 

RESG  = RESG+WG(J) *FSUM 

206 

RESK  = RESK+WGK(JTW) *FSUM 

207 

RESABS  = RESABS+WGK (JTW)*(DABS( F0AL1 ) +DABS ( F0AL2 ) ) 

208 

10 

CONTINUE 

209 

DO  15  J=l, 15 

210 

JTWM1  = J*2-l 

211 

DABSC  = HLGTH*XGK( JTWM1 ) 

212 

CALL  F (CENTR-DABSC, F0AL1 ) 

213 

CALL  F (CENTR+ DABSC, FOAL 2 ) 

214 

FOl(JTWMl)  = FOAL 1 

215 

F02CJTWM1)  = F0AL2 

216 

FSUM  = FOAL 1+F0AL2 

217 

RESK  : RESK+WGKCJTWMl ) *FSUM 

218 

RESABS  = RESABS+WGK  C JTWM1 )*(DABS(F0AL1 )+DABS(F0AL2) 

219 

15 

CONTINUE 

220 

RESKH  = RESK*0. 5D+00 

221 

RESASC  = WGK (31 ) *DABS ( FC- RESKH ) 

222 

DO  20  J = 1 , 30 

223 

RESASC  = RESASC+WGK ( J ) * ( DABS ( FO 1 ( J ) -RESKH )+DABSCF02 

224 

20 

CONTINUE 

225 

RESULT  r RESK*HLGTH 

226 

RESABS  : RESABS* DHLGTH 

227 

RESASC  : RESASC *DHLGTH 

, 228 

ABSERR  = DABS ( ( RESK -RESG ) *HLGTH ) 

91 


I 


229 

230 

231 

232 

233 

234 


I F ( RESASC . NE . 0 . 0D+00 . AND .ABSERR.NE.O. OD+OO) 

* ABSERR  = RESASC*DMIN1 (0. 1D+01 * ( 0. 2D+03*ABSERR/RESASC ) **  1 . 5D+00) 
IF(RESABS.GT.UFL0U/(0.5D  + 02*EPriACH)  > ABSERR  = DMAX1 

* C (EPMACH*0. 5D+02)*RESABS, ABSERR) 

RETURN 

END 


92 


DOUBLE  PRECISION  FUNCTION  D1MACH( I ) 


- 


1 

2 

3 

4 

5 

6 
7 

e 

9 

10 

n 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 


53 

54 

55 

56 

57 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DOUBLE-PRECISION  MACHINE  CONSTANTS 

D1MACHC  1)  = B**<EMIN-1),  THE  SMALLEST  POSITIUE  MAGNITUDE. 

D1MACH ( 2)  = B**EMAX* ( 1 - B**(-T>),  THE  LARGEST  MAGNITUDE. 

D1MACH ( 3)  = B**(-T),  THE  SMALLEST  RELATIUE  SPACING. 

D1MACH ( 4)  = B** ( 1 -T ) / THE  LARGEST  RELATIUE  SPACING. 

D1MACH ( 5)  = LOG 1 0(B) 

TO  ALTER  THIS  FUNCTION  FOR  A PARTICULAR  ENUIRONMENT » 

THE  DESIRED  SET  OF  DATA  STATEMENTS  SHOULD  BE  ACT I UATED  BY 
REMOUING  THE  C FROM  COLUMN  1. 

WHERE  POSSIBLE,  OCTAL  OR  HEXADECIMAL  CONSTANTS  HAUE  BEEN  USED 
TO  SPECIFY  THE  CONSTANTS  EXACTLY  WHICH  HAS  IN  SOME  CASES 
REQUIRED  THE  USE  OF  EQUIUALENT  INTEGER  ARRAYS. 

INTEGER  SMALL ( 4 ) 

INTEGER  LARGE ( 4 ) 

INTEGER  RIGHTC4) 

INTEGER  DIUER ( 4 ) 

INTEGER  L0G10C4 ) 

DOUBLE  PRECISION  DMACH ( 5 ) 

EQUI UALENCE  C DMACH < 1 ), SMALL < 1 ) ) 

EQUIUALENCE  < DMACH ( 2 ), LARGE ( 1 ) ) 

EQUI UALENCE  ( DMACH ( 3 ) , RI GHT ( 1 ) ) 

EQUIUALENCE  ( DMACH ( 4 ), DIUER ( 1 ) ) 

EQUIUALENCE  ( DMACH ( 5 ) , LOG 1 0 ( 1 ) ) 

MACHINE  CONSTANTS  FOR  THE  BURROUGHS  1700  SYSTEM. 


38 

C 

39 

C 

DATA 

SMALL ( 1 ) 

/ 

ZC00800000 

40 

C 

DATA 

SMALL ( 2 ) 

/ 

2000000000 

/ 

41 

C 

42 

C 

DATA 

LARGE  < 1 ) 

/ 

ZDFFFFFFFF 

/ 

43 

C 

DATA 

LARGE  C 2 ) 

X 

ZFFFFFFFFF 

/ 

44 

C 

45 

C 

DATA 

RIGHT ( 1 ) 

/ 

ZCC5800000 

/ 

46 

C 

DATA 

RIGHT (2 ) 

/ 

ZOOOOOOOOO 

/ 

47 

c 

48 

c 

DATA 

DIUERC 1 ) 

/ 

ZCC6800000 

/ 

49 

c 

DATA 

DIUER(2> 

✓ 

ZOOOOOOOOO 

/ 

50 

c 

51 

c 

DATA 

L0G1 0 ( 1 ) 

/ 

ZD00E730E7 

/ 

52 

c 

DATA 

L0G1 0 ( 2 ) 

/ 

ZC77800DC0 

/ 

MACHINE  CONSTANTS  FOR  THE  BURROUGHS  5700  SYSTEM. 

DATA  SMALL ( 1 ) ✓ 01771000000000000  / 

DATA  SMALL ( 2 ) / 00000000000000000  / 


93 


50 

. C 

59 

C 

DATA 

LARGE ( 1 ) 

/ 

00777777777777777  / 

60 

c 

DATA 

LARGE (2) 

/ 

00007777777777777  / 

61 

c 

62 

c 

DATA 

RIGHT(1 ) 

/ 

01461000000000000  / 

63 

c 

DATA 

RIGHT ( 2 ) 

/ 

00000000000000000  / 

64 

c 

65 

c 

DATA 

DIUER  < 1 ) 

/ 

01^51000000000000  / 

66 

c 

DATA 

DIUERC2) 

/ 

ooooooooooooooooo  / 

67 

c 

60 

c 

DATA 

L0G1 0 ( 1 ) 

/ 

01157163034761674  / 

69 

c 

DATA 

L0G1 0 ( 2 ) 

/ 

00006677466732724  / 

70 

c 

71 

c 

MACHINE  CONSTANTS  FOR  THE  BURROUGHS  6700/7700 

72 

c 

73 

c 

DATA 

SMALL ( 1 ) 

/ 

01771000000000000  / 

74 

c 

DATA 

SMALL ( 2 ) 

/ 

07770000000000000  / 

75 

c 

76 

c 

DATA 

LARGE ( 1 ) 

/ 

00777777777777777  / 

77 

c 

DATA 

LARGE ( 2 ) 

/ 

07777777777777777  / 

70 

c 

79 

c 

DATA 

RIGHT ( 1 ) 

✓ 

01461000000000000  / 

00 

c 

DATA 

R I GHT ( 2 ) 

/ 

OOOOOOOOOOOOOOOOO  / 

01 

c 

02 

c 

DATA 

D I UER ( 1 ) 

/ 

01451000000000000  / 

03 

c 

DATA 

DIUERC2) 

/ 

ooooooooooooooooo  / 

04 

c 

05 

c 

DATA 

L0G1 0(1 ) 

✓ 

01157163034761674  / 

06 

c 

DATA 

L0G1 0 ( 2 ) 

/ 

00006677466732724  / 

87 

c 

80 

c 

MACHINE  CONSTANTS  FOR  THE  CDC  6000/7000  SERIES 

89 

c 

90 

c 

DATA 

SMALL ( 1 ) 

/ 

006 04 0000 00 0000 00 000 B / 

91 

c 

DATA 

SMALL ( 2 ) 

/ 

OOOOOOOOOOOOOOOOOOOOB  / 

92 

c 

93 

c 

DATA 

LARGE ( 1 ) 

/ 

37767777777777777777B  / 

94 

c 

DATA 

LARGE ( 2 ) 

/ 

371 67777777777777777 B / 

95 

c 

96 

c 

DATA 

R I GHT ( 1 ) 

/ 

15604 000000000 0000 00 B / 

97 

c 

DATA 

R I GHT ( 2 ) 

/ 

15 000000000000000000 B / 

98 

c 

99 

c 

DATA 

D I UER ( 1 ) 

/ 

156 1 4 0000 00 000 000 000 B / 

100 

c 

DATA 

D I UER ( 2 ) 

/ 

1 50 1 OOOOOOOOOOOOOOOOB  / 

101 

c 

102 

c 

DATA 

L0G10 ( 1 ) 

/ 

17164642023241 175717B  / 

103 

c 

DATA 

L0G1 0 ( 2 ) 

/ 

16367571421742254654B  / 

104 

c 

105 

c 

MACHINE  CONSTANTS  FOR  THE  CRAY  1 

106 

c 

107 

c 

DATA 

SMALL (1  ) 

/ 

200004 000000 00 OOOOOOOB  / 

108 

c 

DATA 

SMALL (2) 

/ 

OOOOOOOOOOOOOOOOOOOOB  / 

109 

c 

110 

c 

DATA 

LARGE ( 1 ) 

/ 

577777777777777 777777 B / 

111 

c 

DATA 

LARGE ( 2 ) 

/ 

000007 777 77 7 77 7 77 7 77 7 B / 

112 

c 

113 

c 

DATA 

RIGHTC 1 ) 

/ 

37721 4 00 000 000 000 000 OB  / 

,114 

c 

DATA 

R I GHT ( 2 ) 

/ 

OOOOOOOOOOOOOOOOOOOOOB  / 

94 


115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 


138 

139 

140 

141 

142 

143 

144 

145 

146 

147 


155 

156 

157 

158 

159 


DATA  DIUER(l)  / 377224000000000000000B  / 

DATA  DIUER ( 2 ) / OOOOOOOOOOOOOOOOOOOOOB  / 

DATA  LOGlO(l)  / 377774642023241 175717B  / 

DATA  L0G1 0(2)  / 00000757 1 42 1 742254654B  / 

MACHINE  CONSTANTS  TOR  THE  DATA  GENERAL  ECLIPSE  S/200 

NOTE  - IT  MAY  BE  APPROPRIATE  TO  INCLUDE  THE  FOLLOWING  CARD 
STATIC  DMACH ( 5 ) 

DATA  SMALL/20K, 3*0/ , LARGE/7777 7K, 3*177777K/ 

DATA  RIGHT/31420K, 3*0/ , D I UER/32020K, 3*0/ 

DATA  L0G10/40423K/ 42023K, 50237K, 74776K/ 

MACHINE  CONSTANTS  FOR  THE  HARRIS  220 


133 

C 

DATA 

SMALL ( 1 ) , SMALL ( 2 ) 

/ 

'20000000, 

'00000201 

/ 

134 

C 

DATA 

LARGE ( 1 ) , LARGE (2) 

/ 

'37777777, 

'37777577 

/ 

135 

C 

DATA 

R I GHT (1), RIGHTO) 

/ 

'20000000, 

'00000333 

/ 

136 

C 

DATA 

D I UER ( 1 ) , DIUERO) 

/ 

'20000000, 

'00000334 

/ 

137 

C 

DATA 

L0G10  < 1 ) , LOG 1 0(2) 

/ 

'23210115, 

' 10237777 

/ 

MACHINE  CONSTANTS  FOR  THE  HONEYWELL  600/6000  SERIES. 

DATA  SMALL( 1 ) , SMALL ( 2 ) / 0402400000000,  0000000000000  / 

DATA  LARGE( 1 ) , L ARGE ( 2 ) / 0376777777777,  0777777777777  / 

DATA  RIGHT ( 1 ) , RIGHT (2)  / 0604400000000,  OOOOOOOOOOOCO  / 
DATA  DIUER( 1 ) , DIUERO)  / 0606400000000,  0000000000000  / 

DATA  L0G10( 1 ) ,L0G10(2)  / 0776464202324,  0117571775714  / 

MACHINE  CONSTANTS  FOR  THE  HP  2100 


148 

C 

3 WORD  DOUBLE 

PRECISION 

OPTION  WITH 

FTN4 

149 

C 

150 

C 

DATA 

SMALL ( 1 ), 

SMALL ( 2 ) , 

SMALL ( 3 ) 

/ 

40000B, 

0, 

1 

/ 

151 

C 

DATA 

LARGE (1 ), 

LARGE ( 2 ) , 

LARGE ( 3 ) 

/ 

77777B , 

177777B, 

177776B 

/ 

152 

C 

DATA 

R I GHT ( 1 ) , 

RIGHT O ) , 

RIGHTO) 

/ 

40000B, 

0, 

265B 

/ 

153 

C 

DATA 

D I UER ( 1 ) , 

DIUERO), 

DIUERO) 

/ 

40000B , 

0, 

276B 

/ 

154 

c 

DATA 

L0G1 0 ( 1 ) , 

LOG 10(2), 

LOG 10(3) 

/ 

46420B, 

46502B , 

77777B 

/ 

MACHINE  CONSTANTS  FOR  THE  HP  2100 
4 WORD  DOUBLE  PRECISION  OPTION  WITH  F'TN4 


160 

C 

DATA 

SMALL (1  ), 

SMALL ( 2 ) 

/ 

40000B, 

0 

/ 

161 

C 

DATA 

SMALL  C3) , 

SMALL ( 4 ) 

/ 

0, 

1 

/ 

162 

c 

DATA 

LARGE ( 1 ) , 

LARGE  C 2 ) 

/ 

77777B , 

177777B 

/ 

163 

c 

DATA 

LARGE ( 3 ) , 

LARGE ( 4 ) 

/ 

177777B, 

177776B 

/ 

164 

c 

DATA 

R I GHT ( 1 ) , 

RIGHT ( 2 ) 

/ 

40000B , 

0 

/ 

165 

c 

DATA 

RIGHTO)  , 

RIGHT ( 4 ) 

/ 

0, 

225B 

/ 

166 

c 

DATA 

DIUERC 1 ) , 

D IUER ( 2 ) 

/ 

40000B , 

0 

/ 

167 

c 

DATA 

DIUERO)  , 

DIUER ( 4 ) 

/ 

0, 

227B 

/ 

168 

c 

DATA 

L0G1 0(1 ) , 

LOG  10(2) 

/ 

46420B , 

46502B 

/ 

169 

c 

DATA 

LOGIOO)  , 

LOG 10(4) 

/ 

76747B, 

176377B 

/ 

170 

c 

171 

c 

MACH  I 

NE  CONSTANTS  FOR  THE 

IBM  360/370  SERIES 

95 


172 

173 

174 

173 

176 

177 

178 

179 

180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

21  1 

212 

213 

214 

215 

216 

217 

218 

219 

220 

221 

222 

223 

224 

225 

226 

227 

228 


non 


THE  XEROX  SIGMA  5/7/9,  THE  SEL  SYSTEMS  85/86  AND 
THE  PERKIN  ELMER  (INTERDATA)  7/32. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


DATA  SMALL (1 
DATA  LARGE ( 1 
DATA  RIGHTd 
DATA  DIUER ( 1 
DATA  L0G1 0 ( 1 


) , SMALL ( 2 ) 
) , LARGE ( 2 ) 
) ,RIGHT(2) 
) , DIUER (2 ) 
) , LOG 1 0(2) 


/ Z00100000, 
/ Z7FFFFFFF, 
/ Z33 100000 / 
/ Z34100000, 
/ Z4 1 1 3441 3 , 


ZOOOOOOOO  / 
ZFFFFFFFF  / 
ZOOOOOOOO  / 
ZOOOOOOOO  / 
Z509F79FF  / 


MACHINE  CONSTANTS  FOR  THE  PDP-10  ( KA  PROCESSOR). 


DATA  SMALL (1),SMALL(2)  / "033400000000, 
DATA  LARGE(1),LARGE(2)  / "377777777777, 
DATA  RIGHT(1),RI GHT ( 2 ) / "113400000000, 
DATA  DIUERd), DIUER(2)  / "114400000000, 
DATA  L0G10( 1 ) ,L0G10(2)  / "177464202324, 


"000000000000  / 
"344777777777  / 
"000000000000  / 
"000000000000  / 
"144117571776  / 


MACHINE  CONSTANTS  FOR  THE  PDP-10  (KI  PROCESSOR). 


DATA  SMALL (1), SMALL (2)  / "000400000000, 
DATA  LARGE (1),LARGE(2)  / "377777777777, 
DATA  RIGHT( 1 ) , RIGHT(2)  / "103400000000, 
DATA  DIUERd  ),DIUER(2)  / "104400000000, 
DATA  LOGlOd  ) ,L0G10(2)  / "177464202324, 


"000000000000  / 
"377777777777  / 
"000000000000  / 
"000000000000  / 
"476747767461  / 


MACHINE  CONSTANTS  FOR  PDP-11  FORTRAN'S  SUPPORTING 
32-BIT  INTEGERS  (EXPRESSED  IN  INTEGER  AND  OCTAL). 


c 

DATA 

SMALL (1 ) 

, SMALL (2) 

/ 

8388608,  0 

/ 

c 

DATA 

LARGEd) 

, LARGE ( 2 ) 

/ 

2147483647,  -1 

/ 

c 

DATA 

RIGHT ( 1 ) 

, R I GHT ( 2 ) 

/ 

612368384,  0 

/ 

c 

DATA 

DIUERd  ) 

, D I UER ( 2 ) 

/ 

620756992,  0 

/ 

c 

c 

DATA 

LOGlOd  ) 

, L0G1 0 ( 2 ) 

/ 

1067065498,  -2063872008 

/ 

Smt 

c 

DATA 

SMALL (1 ) 

, SMALL ( 2 ) 

/ 

000040000000,  000000000000 

/ 

c 

DATA 

LARGEd  ) 

, LARGE ( 2 ) 

/ 

017777777777,  037777777777 

/ 

c 

DATA 

RIGHT (1 ) 

, RIGHT (2 ) 

/ 

004440000000,  000000000000 

/ 

c 

DATA 

DIUERd  ) 

, D I UER ( 2 ) 

/ 

004500000000,  000000000000 

/ 

c 

r* 

DATA 

LOGlOd  ) 

, L0G1 0 ( 2 ) 

/ 

007746420232,  020476747770 

/ 

c 

MACHINE  CONSTANTS  FOR 

PDP 

-11  FORTRAN'S  SUPPORTING 

c 

c 

16-BIT  INTEGERS  (EXPRESSED  IN  INTEGER  AND  OCTAL). 

V 

C 

DATA 

SMALL ( 1 ) 

, SMALL ( 2 ) 

/ 

128,  0 / 

C 

r 

DATA 

SMALL ( 3 ) 

, SMALL ( 4 ) 

/ 

0,  0 / 

L 

C 

DATA 

LARGEd  ) 

, LARGE ( 2 ) 

/ 

32767,  -1  / 

C 

r 

DATA 

LARGE ( 3 ) 

, LARGE ( 4 ) 

/ 

-1,  -1  / 

V 

C 

DATA 

RIGHTd  ) 

, R I GHT ( 2 ) 

/ 

9344,  0 / 

C 

DATA 

RIGHT (3) 

, R I GHT ( 4 ) 

/ 

0,  0 / 

c 

DATA 

DIUERd  ) 

, DIUER(2) 

/ 

9472,  0 / 

c 

r 

DATA 

DIUER ( 3 ) 

, DIUER(4) 

/ 

0,  0 / 

c 

DATA 

LOGlOd  ) 

, LOG 1 0(2) 

/ 

16282,  8346  / 

c 

DATA 

L0G1 0 ( 3 ) 

, L0G10 ( 4 ) 

/ 

-31493,  -12296  / 

96 


230 

C 

DATA 

SMALL ( 1 ) , SMALL ( 2 ) 

/ 

0000200.  0000000  / 

231 

C 

DATA 

SMALL ( 3 ) » SMALL ( 4 ) 

✓ 

0000000.  0000000  / 

232 

C 

233 

C 

DATA 

LARGE (1 ) » LARGE ( 2 ) 

✓ 

0077777.  0177777  / 

234 

C 

DATA 

LARGE(3),  LARGE ( 4 ) 

/ 

0177777.  0177777  ✓ 

235 

C 

236 

C 

DATA 

RIGHT(1),RIGHT(2) 

/ 

0022200.  0000000  / 

237 

c 

DATA 

RIGHTO)  ,RIGHT(4) 

/ 

0000000.  0000000  / 

238 

c 

239 

c 

DATA 

DIUER d ) , DIUER ( 2 ) 

/ 

0022400.  0000000  / 

240 

c 

DATA 

DIUER ( 3 ) / DIUER ( 4 ) 

/ 

0000000.  ooooooo  / 

241 

c 

242 

c 

DATA 

LOGlOd  ) » L0G1 0 ( 2 ) 

/ 

0037632,  0020232  ✓ 

243 

c 

DATA 

L0G10O)»L0G10(4) 

/ 

0102373,  0147770  / 

244 

c 

245 

c 

MACHINE  CONSTANTS  FOR 

THE 

UNIUAC  1100  SERIES.  FTN  COMPILER 

246 

c 

247 

c 

DATA 

SMALL  < 1 ) * SMALL ( 2 ) 

/ 

0000040000000,  0000000000000 

/ 

248 

c 

DATA 

LARGE ( 1 ) , LARGE ( 2 ) 

/ 

0377777777777,  0777777777777 

/ 

249 

c 

DATA 

RIGHT d ) , RIGHT (2 ) 

/ 

0170540000000,  0000000000000 

/ 

250 

c 

DATA 

DIUERd ) , DIUER ( 2 ) 

/ 

0170640000000,  0000000000000 

/ 

251 

c 

DATA 

L0G10< 1 ) , LOG 10 (2 ) 

/ 

0177746420232,  0411757177572 

/ 

252 

c 

253 

c 

254 

c 

MACHINE  CONSTANTS  FOR 

THE 

UNIUAC  1100  SERIES.  FOR  COMPILER 

255 

c 

256 

c 

DATA 

SMALL( 1 ) ,SMALL(2) 

✓ 

0000040000000,  0000000000000 

/ 

257 

c 

DATA 

LARGEd  ) , LARGE  ( 2 ) 

/ 

0377777777777,  0777777777777 

/ 

258 

c 

DATA 

RIGHT  d ) .RIGHTO) 

/ 

0170540000000,  0000000000000 

/ 

259 

c 

DATA 

DIUERd  ) , DIUERO) 

/ 

0170640000000,  0000000000000 

/ 

260 

c 

DATA 

LOGlOd  ) ,L0G10(2) 

/ 

0177746420232,  0411757177572 

/ 

261 

c 

262 

c 

MACHINE  CONSTANTS  FOR 

UAX 

-11/780 

263 

c 

(EXPRESSED  IN  INTEGER 

AND 

OCTAL) . 

264 

c 

265 

c 

DATA 

SMALL(l).SMALLO) 

/ 

128,  0 / 

266 

c 

DATA 

LARGEd  ) » LARGE  ( 2 ) 

/ 

-32769,  -1  / 

267 

c 

DATA 

RIGHT ( 1 ) » R I GHT ( 2 ) 

/ 

9344,  0 ✓ 

268 

c 

DATA 

DIUER  ( 1 ) , DIUERO) 

/ 

9472,  0 / 

269 

c 

DATA 

LOGlOd  ) .LOGIOO) 

/ 

546979738,  -805665541  / 

270 

c 

271 

c 

DATA 

SMALL ( 1 ) . SMALL ( 2 ) 

/ 

000000000200,  000000000000  / 

272 

c 

DATA 

LARGEd  ) » LARGE  ( 2 ) 

/ 

037777677777,  037777777777  / 

273 

c 

DATA 

RIGHTd  ), RIGHTO) 

/ 

000000022200,  000000000000  / 

274 

c 

DATA 

DIUERd  ) , DIUERO) 

/ 

000000022400,  000000000000  / 

275 

c 

DATA 

LOGlOd  ). LOGIOO) 

/ 

004046437632,  031776502373  / 

276 

c 

277 

c 

IF  (I  .LT.  1 .OR.  I 

.GT 

. 5) 

278 

c 

1 CALL  SETERR( 24HD1MACH 

- I OUT  OF  BOUNDS, 24, 1, 2) 

280 

281 

282 

283 


D1MACH  = DMACH(I) 
RETURN 

END 
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STITL  PWAFFI  — Initialize  form-factor  table  for  C-12 
SUBROUTINE  PWAFFI 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C The  size  of  the  arrays  in  the  common  PWAFF 

PARAMETER  (MNPT=200) 

C The  number  of  points  in  the  interpolation  table 

PARAMETER  <NTBM1=179) 

C The  number  of  points  + 1 

PARAMETER  (NTB=NTBM1+1 ) 

DOUBLE  PRECISION  XC 1 ( 1 00 ) , XC2 ( 79 ) , FC 1 C 1 00 ) , FC2 < 79 ) , X ( MNPT ) 
DOUBLE  PRECISION  F ( NTBM1 ), DF < NTBM1 ), A ( MNPT ), B < MNPT ), C ( MNPT ) 
DOUBLE  PRECISION  D ( MNPT ), AP ( MNPT ), BP ( MNPT ), CP ( MNPT ), DP C MNPT ) 
DOUBLE  PRECISION  XC ( NTBM 1 ) , FC C NTBM 1 ) 

COMMON  /CONST / EM, EM2 , AM , HB ARC , ALPH, P IE , ALPHA , RTD , DTR 
COMMON  /PWAFF/  NPT, X, A, B, C, D, AP, BP, CP, DP 
COMMON  /CORREC/IWTILD, WTIL1 , WTIL2, I SOFT, S0FT1 
EGUIUALENCE  ( FC ( 1 ) , FC1 ( 1 ) ) , ( FC ( 101 ) , FC2 ( 1 ) ) 

EQUIUALENCE  ( XC ( 1 ) , XC 1 ( 1 ) ) , ( XC < 1 01 ) , XC2( 1 ) ) 

Table  of  q and  corresponding  form  factor  from  DWBA  code.  In 
principle  this  table  could  be  read  in  rather  than  defined  in 
data  statements. 

DATA  XC1/5. 3068D-02, 1 . 06 1 3D-01 , 1 . 5918D-01, 2. 1 222D-0 1 , 2 . 6523D-01 
+ 3. 1 822D-01 , 3. 71 18D-01, 4. 2410D-01, 4. 7697D-01 , 5. 2980D-01 , 

+ 5. 8258D-01, 6.3530D-01,6. 8795D-01 , 7 . 4054D-01 , 7.9305D-01 , 

+ 8 . 4549D-0 1 »8.97B4D-01,9.5010D-01, 1.0023D+00, 1.0543D+00, 

+ 1 . 1 063D+00,  1 . 1582D  + 00, 1 . 2099D+00, 1 . 2615D  + 00,  1 . 3130D  + 00, 

+ 1 . 3644D+00, 1 . 4156D+00, 1 . 4667D+00, 1 . 5177D+00, 1 . 5685D+00, 

+ 1 . 6191 D + 00, 1 . 6696D+00,  1 . 7199D  + 00,  1 . 7701 D+00,  1 . 8201 D+00, 

+ 1 . 8699D  + 00,  1 . 9195D+00, 1 . 9690D+00, 2. 0182D  + 00, 2. 0673D+00, 

+ 2.1 162D  + 00, 2. 1648D+00, 2 . 2 1 33D  + 00 , 2 . 26 1 5D+ 00 , 2 . 3096D  + 00 , 

+ 2. 3574D  + 00, 2. 4050D  + 00, 2. 4524D+00, 2. 4 995D+00,  2 . 54  64 D + 00, 

+ 2. 5931 D+00, 2. 6395D  + 00, 2. 6856D  + 00, 2. 7316D+00,  2. 7772D  + 00, 

+ 2. 8226D+00, 2. 8678 D + 00, 2. 9127D  + 00, 2. 9573D+00,  3. 0017D  + 00, 

+ 3. 0457 D + 00, 3. 0895D  + 00, 3. 1331 D + 00, 3. 1763D  + 00,  3. 2193D  + 00, 

+ 3. 2619D  + 00, 3 . 3043D  + 00, 3. 34  64D  + 00,  3.  3881 D + 00,  3. 4296D  + 00, 

+ 3. 4708D  + 00, 3. 51 16D+00, 3. 5522D  + 00, 3. 5924D  + 00,  3. 6324 D + 00, 

+ 3. 6720D  + 00, 3. 71 13D+00, 3 . 7502D  + 00, 3. 7889D  + 00,  3. 8272D  + 00, 

+ 3. 8652D  + 00, 3. 9029D  + 00, 3. 9402D  + 00, 3. 9772D+00,  4 . 0138D  + 00, 

+ 4. 0501 D+00, 4. 0861 D+00, 4. 1218D+00,4. 1570D+00,4. 1920D+00, 

+ 4. 2266D+00, 4 . 2608D+00, 4 . 2947D+00, 4. 3282D+00, 4. 36 14 D+00, 

+ 4 . 3942D+00, 4 . 4267D+00, 4 . 4588D+00, 4. 4906D+00, 4 . 5219D+00/ 

DATA  XC2/4. 5529D+00,4. 5836D400, 4. 6139D+00, 4 . 6438D+00, 4 . 6733D+00 
+ 4. 7025D+00, 4 . 7313D+00, 4 . 7598D+00, 4 . 7878D+00, 4 . 8155D+00, 

+ 4. 8428D+00, 4 . 8698D+00, 4 . B963D+00, 4. 9225D+00, 4. 9483D+00, 

+ 4. 9737D+00, 4 . 9987D  + 00, 5 . 0234D+00, 5. 0477D+00,  5 . 0715D  + 00, 

+ 5. 0950D+00, 5. 1182D+00,5. 1409D+00,5. 1632D+00,5. 1852D+00, 

+ 5. 2068D+00, 5. 2280D+00, 5. 24  88D  + 00, 5. 2692D  + 00,  5. 2892D  + 00, 

+ 5. 3088D+00, 5. 3280D+00, 5. 3469D+00, 5 . 3653D+00, 5. 3834D+00, 

+ 5. 4011 D+00, 5 . 4183D+00, 5. 4352D+00, 5. 4517D+00, 5. 4678D+00, 

+ 5. 4935D+00, 5 . 4988D+00, 5 . 51 37D+00, 5. 5283D+00, 5. 5424 D+00, 

+ 5. 5561 D + 00, 5. 5695 D + 00, 5. 5824D  + 00, 5. 5949D+00,  5. 6071 D + 00, 

+ 5. 61B8D  + 00, 5. 6302D  + 00, 5. 641 1 D + 00, 5. 6517D+00,  5. 6618D  + 00, 

+ 5. 6716D+00, 5. 6810D+00, 5. 6899D+00, 5. 6985D+00, 5 . 7067D+00, 

+ 5. 7 145D+00, 5.72 18D+00, 5 . 7288D+00, 5. 7354 D+00, 5. 7416D+00, 

+ 5. 7474 D+00, 5 . 7527D  + 00 , 5 . 7577D  + 00,  5. 7623D  + 00,  5 . 7665D+00, 

+ 5. 7703D+00, 5. 7737D  + 00,5. 7767D  + 00,  5. 7793D+00,  5. 7815D  + 00, 
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59 
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+ 5. 7833D+00, 5. 7847D+00, 5 . 7857D+00 , 5 . 7862D+00/ 

DATA  FC1/1 . 0000D  + 00 > 9 . 8087D-01, 9. 5239D-01, 9. 1556D-01 / 8. 697  5D-01 , 
+ 8. 1646D-01, 7. 5736D-01 , 6. 9420D-01, 6. 2873D-01, 5. 6261 D-01, 

+ 4. 9734D-01, 4.3427D-01, 3. 7450D-01 , 3. 1890D-01 , 2 . 6807 D-01  , 

+ 2 . 2240D-0 1 , 1 . 8203D-01 , 1 . 4692D-01, 1 . 1688D-01, 9. 1572D-02, 

+ 7. 0602D-02, 5. 3507D-02,3. 9805D-02, 2. 9015D-02, 2. 0675D-02, 

+ 1. 4357D-02,  9. 6756D-03, 6. 2920D-03, 3 . 916 1 D-03, 2 . 3048D-03, 

+ 1 . 2588D-03, 6. 1895D-04, 2 . 6065D-04, 8. 91 12D-05, 3 . 4238D-05, 

+ 4 . 5899D-05, 8.9708D-05, 1 . 4332D-04,  1 . 9333D-04,  2. 3265D-04, 

+ 2. 5852D-04, 2. 7093D-04, 2. 71 43D-04 , 2 . 6234D-04,  2 . 4618D-04 , 

+ 2. 2533D-04, 2. 0186D-04, 1 . 7744D-04, 1 . 5334 D-04, 1 . 3047 D-04, 

+ 1 . 0941 D-04, 9. 0497D-05,7. 3864D-05 , 5 . 951 2D-05 , 4. 7338D-05, 

+ 3. 7177D-05,  2.8825D-05/2. 2064D-05, 1 . 6672D-05,  1.2435D-05, 

+ 9. 1554D-06, 6.6544D-06, 4 . 7753D-06, 3 . 3842D-06, 2 . 3691 D-06 , 

+ 1. 6387D-06,  1 . 1202D-06»7. 56 6 2D -07 , 5. 0471D-07, 3. 3205D-07  » 

+ 2. 1493D-07,  1 . 3630D-07, 8. 4 1 39D-08 , 5 . 0045D-08 , 2 . 825 1 D-08 , 

+ 1 . 4815D— 08, 7.0487D-09, 3. 1135D-09, 1 . 7327D-09, 2 . 0021 D-09 , 

+ 3. 2652D-09,  5.0351D-09, 6. 9478D-09, 8 . 7363D-09, 1 . 0215D-08, 

+ 1 . 1272D-08,  1 . 1 855D-08 » 1 . 1969D-08,  1 . 1660D-08,  1 . 1001D-08, 

+ 1 . 0085D-08,  9. 0060D-09, 7 . 8529D-09, 6. 7016D-09, 5 . 6109D-09, 

+ 4. 6210D-09, 3. 7547D-09, 3. 0203D-09, 2. 4145D-09, 1 . 9268D-09/ 

DATA  FC2/1 . 5422D-09, 1 . 2443D-09, 1 . 0168D-09,  8 . 4 509D-1 0 , 7 . 1668D-10, 
+ 6. 2124D-10, 5.5047D-10/4. 9786D-10, 4. 5831D-10,  4. 2782D-10, 

+ 4. 033 ID- 10, 3. 8236D-10, 3. 6319D-10, 3. 4448D-10, 3. 2538D-10, 

+ 3. 054 ID- 10, 2. 8440D-10, 2. 6246D-10, 2. 3986D-10, 2. 1698D-10, 

+ 1 . 9427D-10, 1 . 7217D— 10, 1.5108D-10, 1.3133D-10, 1. 1315D-10, 

+ 9.6708D-11, 8. 206 ID- 11, 6. 9205D-11, 5.8070D-11, 4.B544D-11, 

+ 4. 0484D-11, 3. 3732D-1 1 ,2. 8123D-11 , 2. 3498D-1 1,  1 . 9707D-1 1 , 

+ 1.6614D-11, 1 .4098D-11, 1 .2055D-11, 1 . 0398D- 1 1 , 9 . 05 1 4 D- 1 2 , 

+ 7. 9547D-12, 7.0577D-12, 6. 31 95D- 1 2 , 5 . 707 1 D- 1 2 , 5 . 1943D-12, 

+ 4. 760 ID- 12, 4. 3877 D- 12, 4 . 0639D-12, 3. 7782D-12, 3. 5222D-12, 

+ 3. 2896D-12, 3. 0756D-12, 2. 8765D-12, 2. 6900D-12, 2. 5143D-12, 

+ 2. 3433D-12, 2. 1 91 1 D-12, 2 . 0422D-12, 1 . 9009D-12, 1 . 7673D-12, 

+ 1 . 641 2D- 12, 1 .5229D-12, 1 . 4126D-12, 1 . 3105D-12, 1 . 2164D-12, 

+ 1. 1 302D-1 2, 1.051 3D- 12, 9. 792 ID- 13, 9. 1 362D- 1 3, 8 . 5462D- 1 3 , 

+ 8. 026 ID- 13, 7 .5814 D- 13, 7 . 21 58D- 1 3 , 6 . 9264D- 1 3 , 6 . 6982D- 1 3 , 

+ 6. 4928D-13, 6. 2152D-13, 5. 5475D-13, 3. 2798D-13/ 

Convert  table  from  form  factor  versus  q to  form  factor 
times  2**2  versus  q**2 
DO  10  1=1, NTBM1 

X ( I ) =XC( I ) * *2 
F ( I )=FC( I ) *36 . DO 

IF  ( IWTILD . NE . 0 ) F C I ) =F < I ) * ( WT I L 1 +WT IL2*DL0G ( X ( I ) /EM2 ) ) 

10  CONTINUE 

C 

C Find  the  cubic  coefficients  for  spline  interpolation  of 

C the  table.  A good  reference  for  spline  interpolation  is 

C Elementary  numerical  analysis  by  S.D.  Conte  and  Carl  deBoor. 

A ( 1 ) =36 . DO 
B( 1 ) =0. DO 
CC 1 ) =0. DO 
D( 1 ) =0. DO 
X(NTB) =1 . D74 
A ( NTB ) = 0 . DO 
BCNTB)=0. DO 
C(NTB) =0. DO 
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DCNTB) =0. DO 

C Calculate  the  derivatives  of  the  form  factor  ( DF ) according 

C to  the  spline  method  at  each  of  the  tabulated  points. 

P=0. DO 
QrO. DO 

H2r 1 .D0/(X(2)-X( 1) ) 

DO  908  I z 2 / NTBM 1 
H1=H2 
H12=H1*H1 

H2= 1 . DO/(X( 1+1 ) -X ( I ) ) 

H22zH2*H2 

Rz2. D0*(H1+H2)-P 

B ( I ) :H2/R 

DF< I ) z-3. DO* ( F ( 1-1 )*H12+F( I ) * ( H22-H12 ) -F ( I +1 )*H22) 
DF(I)=(DF(I)-Q)/R 
PrB( I )*H2 
0=  DF ( I )*H2 

908  CONTINUE 

C Now  back  substitution 

PrO. DO 

DO  909  I r 1 / NTBM1 
JrNTBMl+1-I 
DFC J)zDF( J)-P*B( J) 

PrDFC J) 

909  CONTINUE 

C Using  the  spline  approximation  values  for  the  derivatives 

C at  the  points  in  the  table/  calculate  the  coefficients 

C for  a cubic  which  exactly  fits  these  derivatives  and  the 

C values  of  the  function  at  the  end  points  of  each  interval. 

DO  20  I z2/  NTBM1 

Hr  1 . DO/(X( I )-X( 1-1 ) ) 

H2r H*H 

AlrDF (I)-2. DO*F ( I ) *H 
Cl r F ( I )-Al*X( I ) 

A2zDF ( 1-1 ) +2  . DO*F ( I-i )*H 
C2=F( 1-1 ) -A2*X ( 1-1 ) 

D(I)r(Al+A2)*H2 

C( I >r (Cl -2. D0*A1*X(I-1 )+C2-2. D0*A2*X< I ) )*H2 
B ( I ) z ( A 1 *X ( I - 1 ) **2-2 . D0*C1 *X( 1-1 ) +A2*X( I )**2- 
+ 2. D0*C2*X( I ) )*H2 

A(I)r(Cl*X(I-l) **2+C2*X( I )**2)*H2 
20  CONTINUE 

C For  this  model/  we  assume  that  Wl  (the  transverse  form  factor) 

C is  0. 

DO  30  Irl,NTB 
AP( I ) rO. DO 
BP( I ) rO. DO 
CP( I ) rO. DO 
DPC I ) rO. DO 
30  CONTINUE 

NPTzNTB 
RETURN 
END 
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STITL  F0RM1  — Transverse  form  factor  by  interpolation 
FUNCTION  F0RMKQ2) 

C Size  of  arrays  in  PWAFF 

PARAMETER  (MNPT=200) 

IMPLICIT  DOUBLE  PRECISION  (A-H/O-Z) 

DOUBLE  PRECISION  X ( MNPT) / A ( MNPT ) / B ( MNPT ) » C ( MNPT )<D( MNPT ) 
DOUBLE  PRECISION  AP ( MNPT ),  BP (MNPT ), CP < MNPT ), DP ( MNPT ) 
COMMON  /PWAFF ✓ NPT/X»A»B#C»D»AP»BP#  CP»  DP 
C CALL  STARTT(l) 

0Q=-Q2 

DO  10  1=1/ NPT— 1 

IF  (QQ.LT.X(I))  GO  TO  20 
10  CONTINUE 

I =NPT 

20  F0RM1 = ( ( DP ( I ) *QQ+CP ( I ) )*QQ+BP( I ) )*QQ+AP( I ) 

C CALL  STOPT(l) 

RETURN 

END 
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*TITL  F0Rr^2  — Longitudinal  form  factor  by  interpolation 
FUNCTION  F0RM2C02) 

C Size  of  arrays  in  PWAFF 

PARAMETER  (MNPT=200) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

. DOUBLE  PRECISION  X ( MNPT ) , A ( MNPT ) , B ( MNPT ) , C ( MNPT ) , D < MNPT ) 
DOUBLE  PRECISION  AP ( MNPT ), BP C MNPT ), CP ( MNPT ), DP ( MNPT ) 
COMMON  /"PWAFF  / NPT*X*A,B/C»D*  AP  / BP  / CP  / DP 
C CALL  STARTT(l) 

QQz-QZ 

DO  10  1=1*  NPT  — 1 

IF  ( QQ . LT . X ( I ) ) GO  TO  20 
10  CONTINUE 

I =NPT 

20  F0RM2= ( (DC I)*QQ  + C< I ) >*QQ  + B<  X ) >#QQ  + A( I ) 

C CALL  STOPT(l) 

RETURN 

END 


Manuscrit  raqu  le  7 novembre  1983 
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